NB. Usage: load 'extract_pol.ijs', load 'Set_LMN.ijs', then Set_LMN''
NB. Z=. ('marcs_package';'p3000_g+4.0',tag) ExtractJ 5200 5100 5000
NB. or, with 'lams' a list of wavelengths, e.g., here are 20 lambdas:
NB. lams=: 1e3*8 7.5 7 6.5 6 5.5 5.2 5 4.8 4.6 4.4 4.2 4 3.8 3.65 3.55 3.4 3.2 3 2.8
NB. Z=. ('marcs_package';'sun.opa') Extract lams
NB. The output is: 'mu rad deg lamc Flxs Ia Qa pol'=. Z  
NB. The 1st 3 are the angles: mu & theta(mu) in radians & in degrees
NB. lamc=wavelengths & Flxs are the monochromatic (astrophysical) fluxes
NB. Ia Qa pol are the intensities, Stokes Q, and polarization
NB. where Ia = Ia(lam,mu): $Ia -> 3 159   (159 mu values) 
NB. Try: 'pensize 3' plot mu;pol
require 'files'
NB. !!! note: must run Set_LMN'' first. !!!
 
Extract=: dyad define
nl=. # nus=. 3e18% lams=: y
'fd fn'=. x   NB. fd=file directory, fn=file name
dir=: '/your_path/',fd,'/'
Z=. mread dir,fn       NB. read in the .opa file
AA=. Z do_file lams    NB. <== extract data from file format           
load'/your_path/splines.ijs'
load'/your_path/IQMsi.ijs'
load'/your_path/IQs.ijs'
load'/your_path/s_and_p.ijs'
frame_s=: < (#ttt)$0   NB. ttt is tau grid defined (ttt=: ) in "Set_LMN"
frame_p=: < (#ttt)$0
mua=: 0 0.0001 0.0003 0.001 0.0015 0.002 0.003 0.005 0.0075 0.01 0.013 0.016
theta=. acos mua=: mua, (4}. _101}. int01 200),(50}. int01 100)
E0=. IQMsi ttt;mua  
frame=. < (3,(#mua))$mua, theta, (180*theta%1p1)
ind=. _1
while. (ind=. >:ind)<nl do.        NB. loop over nl wavelengths
'tau T ops abs sca'=. ind{ AA      NB. pull off data for given lam
NB. "abs" and "sca" are scaled: true kappa=abs*ops, sigma=sca*ops
               NB. We add tau=0 point, extrapolating all variables:
T=. T,~ (tau;T) splint 0
ops=. ops,~ (tau;ops) splint 0
abs=. abs,~ (tau;abs) splint 0
sca=. sca,~ (tau;sca) splint 0
tau=. 0, tau
rat=. abs+sca                      NB. ratio of kappa(nu) to standard kappa
taunu=. tau spintg rat             NB. integrate k_nu/k_std to get tau(nu)
nu=. ind{nus                       NB. current frequency
B0=. (T0=. 40{T) Bnu nu            NB. evaluate B_nu at T0 (standard tau=0.566)
B=. B0%~ (nu Bnu~ T)               NB. evaluate B_nu(T) and normalize to B_nu(T0)
BB=. (taunu;B) splint ttt
lam=. abs%rat                      NB. lambda = kappa_abs/(kappa_abs+sigma)
llam=. (taunu;lam) splint ttt
NB. matrix transforms L0,M0,N0 defined externally in "Set_LMN".
's p'=: (L0;M0;N0) s_and_p (llam;BB)  NB. call to solve for s & p  
s=: B0*s                           NB. Undo normalization
p=: B0*p
frame_s=: frame_s,< >s
frame_p=: frame_p,< >p
NB. 'Ia Qa pol'=: (E0;mua) IQs s;p NB. use E-transform to get I&Q from s&p
frame=. frame, < >(E0;mua) IQs s;p
end.
frame_s=: >}. frame_s
frame_p=: >}. frame_p
Flxs=: +/PHI2* |:frame_s
Flxs=: Flxs+ +/PHI4* |:frame_p
Z=. >frame
'mu rad deg'=. 0{Z
'Ia Qa pol'=. |:"2|:|:"2 z1=. }.Z
mu;rad;deg;lams;Flxs;Ia;Qa;pol
)

NB. extract data from MARCS formated .opa file
do_file=: dyad define  
ndp=. 1{0{x      NB. # of depth points
swave=. 2{0{x    NB. standard wavelength
nwave=. 0{1{x    NB. # of wavelengths 
Z=. 2}.x         NB. trim off 1st 2 rows
nwl=. >. nwave % 10   NB. # of wavelength lines
lam=. (i. nwave){ ,(i. nwl){Z
Z=. nwl }. Z          NB. trim off wavelengths
nedp=. >: nwave%3     NB. # rows in each depth block
ndl=. ndp*nedp        NB. number of data rows
Z=. (i. ndl){ Z       NB. data block (20048 x 10)
Z=. (ndp,nedp,10)$,Z  NB. (56 x 358 x 10) array
H=. _2}."1 ({."2 Z)   NB. take 1st row of block & trim last 2 columns 
'rad tau0 T Pe Pg rho xi ops'=. |:H
A=. >tau0;T;ops
B=. _4}."1 (}."2 Z)
B=. (ndp,(<:nedp),3,2)$ ,B
abs=. ,"2 (0{"_3 B)   NB. extract (56 x 1071) absorption array
sca=. ,"2 (1{"_3 B)   NB. scattering array
nww=. #ww=. y         NB. ww are the target wavelengths 
Qa=. (lam;|:abs) lintrp ww  NB. interpolate abs at target wavelength
Qs=. (lam;|:sca) lintrp ww  NB. interpolate sca at target wavelength
A=. (nww,($A))$,A     NB. one set of (tau0,T,Pg,ops) for each ww
(A,"_1 Qa),"_1 Qs     NB. append abs(lam) & sca(lam) to each set
)

NB. linear interpolation: (x;y) lintrp x0 --> y(x0)
NB. x & y boxed, x0 may be a vector, extrapolate off ends.
lintrp=: dyad define
'xx yy'=. x
g0=. (}: yy)% d=. DEL xx
g1=. (}. yy)% d
u=. (g0* }.xx) - g1* }:xx
n=. xx IdX y
(n{u)+ y*n{g1-g0 
)
IdX=: ] I.~ [: }: [: }. [

mread=: (0&".)@('m'&fread)  NB. read in data matrix
NB. Planck function B_nu(T) (not pi*B_nu):
Bnu=: 1.47449934e_47&*@(^&3@])% -&1@(^@(4.799243e_11&*@(] % [)))

NB. last part of name for MARCS models:
tag=: '_m0.0_t01_st_z+0.00_a+0.00_c+0.00_n+0.00_o+0.00_r+0.00_s+0.00.opa'
NB. last part of name for models in "marcs_z-1":
tag3=: '_m0.0_t01_st_z-1.00_a+0.40_c+0.00_n+0.00_o+0.40_r+0.00_s+0.00.opa'
NB. last part of name for models in "marcs_package":
tag4=: '_m0.0_t01_st_z+0.00_a+0.00_c+0.00_n+0.00_o+0.00_r+0.00_s+0.00.opa'
NB. last part of name for models in "marcs_z-2":
tag5=: '_m0.0_t01_st_z-2.00_a+0.40_c+0.00_n+0.00_o+0.40_r+0.00_s+0.00.opa'
NB. last part of name for models in "marcs_hot":
tag6=: '_m0.0_t02_st_z+0.00_a+0.00_c+0.00_n+0.00_o+0.00_r+0.00_s+0.00.opa'