NB. Polarization from HOT stellar atmospheres.
NB. Usage: load 'extract_hot.ijs', load 'Set_LMN.ijs'
NB. run  Set_LMN'' once, then do, e.g., 
NB. Z=. 't380g350' Extract 4000  (for T=35,000, log g=3.5)
NB. The output is: 'mu theta(rad) theta(deg)'=. 0{Z
NB. 'Ia Qa pol'=. 1{Z  (1 is wavelength index)
NB. where Ia = Ia(mu) for (#mua) values of mu
NB. dimensions of "Z" are ((#lams),3,(#mua))                              

require 'files'
NB. note: must run Set_LMN'' first.
Extract=: dyad define
nl=. # nus=. 3e18% lams=: y
dir=. '/your_path/STERNE_HOT.d/'
out_file=. dir,'temp0'
in_file=. dir,x,'p00.summary'
2!:1 'awk ''NR >7'' <',in_file,' >',out_file
NB. 2!:1 ... : execute awk command to strip off 1st 7 lines of file
opN=. mread out_file     NB. read in atmosphere table          
opN=. (nl,50,13)$ ,opN   NB. reformat to separate wavelength groups
print 'Wavelengths: ';lams     
load'/juno/jph/J6/SPLINES/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)
I0=. IQMsi ttt;mua  
frame=. < (3,(#mua))$mua, theta, (180*theta%1p1)
ind=. _1
while. (ind=. >:ind)<nl do.
'n ltau T delT lPg lPe ops sca rho dst j1 j2 j3'=. |: ind{ opN
tau=. 10^ltau
abs=. ops - sig                    NB. the "kap" column is the total opacity
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)%ops                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%(abs+sca)                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 via  LAPACK 
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
>frame
)

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

NB. Here we assume the directory "STERNE_HOT.d" contains the model 
NB. atmosphere files with names like "t380g350p00.summary" from the  
NB. website "http://star.arm.ac.uk/~csj/models/Grid.html"
NB. (these model listings give kappa for a wavelength of 4000A only)