NB. Usage: load 'Extract_Kurucz.ijs' This is to extract the polarization
NB. and limb darkening from Kurucz atmospheres processed through the
NB. synspec51 program, modified to write fort.96 (see "Stellar Polarization"
NB. link on the main page above "J and Astronomy").
NB. The output is: 'mu rad deg lam Flx Ia Qa pol'=. Z =. Extract'' 
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) 
require 'files'
load'/your_path/Set_LMN.ijs'
Extract=: monad define
Set_LMN''
'tauR Irho T P Ne k_R'=. read8''
'mass T Ne rho'=. read6''
z=: 0,+/\ dz=: (DEL Irho)%(BIN rho)
lam=. 1{"1 g=. read96'' 
print ({.{.lam),({.{:lam)
Achi=. 2{"1 g                 NB. this is rho*(kappa+sigma) 
Asca=. 5{"1 g                 NB. this is rho*sigma
erase 'g'
tnu=. |: 0, +/\dz* |: BIN"1 Achi
nl=. # nus=. 3e18% lams=: {."1 lam
load'/your_path/SPLINES/splines.ijs'
load'/your_path/IQMsi.ijs'
load'/your_path/IQs.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)
NB. mua=: 0 0.005 0.01 0.025 0.05 0.075 0.1 0.15 0.2 0.25 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
NB. theta=. acos mua
I0=. IQMsi ttt;mua  
NB. matrix transforms L0,M0,N0 defined externally in "Set_LMN".
c=. (3%8)*M0
d=. (3%8)*N0
QQ=. (L0,"1 M0%3),"2 c,"1 d   NB. partial coefficient matrix QQ=[a,b/c,d]
frame=. < (3,(#mua))$mua, theta, (180*theta%1p1)
ind=. _1
while. (ind=. >:ind)<nl do.
  chi=. ind{Achi                    NB. this is rho*(kappa+sigma)
  sca=. ind{Asca                    NB. this is rho*sigma
  tau=. ind{tnu
 nu=. ind{nus                       NB. current frequency
 B0=. (T0=. 52{T) Bnu nu            NB. evaluate B_nu at T0 (standard tau=0.562)
 B=. B0%~ (nu Bnu~ T)               NB. evaluate B_nu(T) and normalize to B_nu(T0)
 BB=. (tau;B) splint ttt
 lambda=. (chi-sca)%chi             NB. lambda = kappa_abs/(chi=kappa_abs+sigma)
 llam=. (tau;lambda) splint ttt
 's p'=: QQ s_p (llam;BB)       NB. call to solve for s & p (direct inversion)
 s=: B0*s                           NB. Undo normalization
 p=: B0*p
 frame_s=: frame_s,< >s
 frame_p=: frame_p,< >p
 NB. 'Ia Qa pol'=: (I0;mua) IQs s;p
 frame=. frame, < >(I0;mua) IQs s;p
end.
frame_s=: >}. frame_s
frame_p=: >}. frame_p
Flxs=: +/PHI2* |:frame_s
Flxs=: Flxs+ +/PHI4* |:frame_p
frame=. >frame
'mu rad deg'=. 0{frame
frame=. }. frame
'Ia Qa pol'=. |:"2|:|:"2 frame
mu;rad;deg;lams;Flxs;Ia;Qa;pol
)

s_p=: dyad define      NB. solve directly by %.
'lam B'=. y            NB. Plank term B(tau) and lambda(tau)
idn=. =@i. +:n=: #B    NB. idn is 2n x 2n identity matrix
oml=: ,~ (1-lam)       NB. (1-lam_i),(1-lam_i)
Q=. idn - x*oml        NB. the grand matrix
rhs=. (lam*(n$B)),n#0  NB. rhs is the right hand side of the system
sp=. rhs%. Q           NB. solve linear system: rhs %. (coeff matrix)
's p'=. (2,n)$ sp      NB. 1st n elements are "s"; 2nd n are "p"
)

read96=: monad define       NB. read data from synmod output fort.96
vv=. mread'fort.96'         NB. (synmod is modified synspec51 program) 
ND=. >./ idx=. 0{"1 vv      NB. number of depth points = ND
hunk=. ND* BS=. (555{. idx)I. 2    NB. block size
nblock=. <. hunk%~ nline=. {. $vv  
nover=. nline- hunk*nblock  NB. number of lines over last whole block
vv=. (-nover)}. vv          NB. trim off final partial block
vv=. (nblock,ND,BS,6)$,vv   NB. groups --> dimensions
vx=: 2{."2 vv               NB. get 1st two rows of each depth group
vx=: |: ,"2 |:"3 |:"2 |:vx  NB. merge into continuous wavelength file
vx=: 2%~ _2 +/\ vx          NB. average continuum set ends
vv=. 2}."2 vv               NB. trim 1st two rows of each depth group
|: ,"2 |:"3 |:"2 |: vv      NB. merge into continuous wavelength file
)

read8=: monad define        NB. read structure from fort.8 (input Kurucz model)
fort8=. 22}. mread 'fort.8'
nd=. 2{ {. fort8            NB. nd=numb depth pts, usually 72
fort8=.  (<:nd){. (2}. fort8)
'Irho T P Ne k_R'=. |: 5{."1 fort8
tauR=. }. |. 100*10^(-0.125*i. nd)  NB. 8 pts/decade, deepest @ 100
>tauR;Irho;T;P;Ne;k_R
)

read6=: monad define       NB. synspec output from Kurucz model input
f6=. 'm' fread 'fort.6'
NB. table starts 3 lines below line including the word "MASS"
idx=. {. 3+ <.({:$f6)%~ I. 'MASS' E. ,f6  
f6=. 0&". (D2e f6)
ND=. 71              NB. number of depth points
f6=. (idx+i. ND){f6
|: (1 2 3 4){"1 f6
)

D2e=: 3 : '''e''(I. ''D''=y) } y' "0  NB. change Fortran D's to e's
mread=: (0&".)@('m'&fread)  NB. read in a data matrix
Bnu=: 1.47449934e_47&*@(^&3@])% -&1@(^@(4.799243e_11&*@(] % [))) NB. Planck fn