NB. Polarization across line in a slowly rotating (spherical) 
NB. star using MARCS model atmosphere tables. Lambda is the  
NB. wavelength in A. The inclination is fixed at 90 deg. 
NB. (Other inclination angles give identical results except 
NB. the rotational velocity Doppler shift is scaled by sin i.)
NB. This uses a Milne-Eddington line profile to evaluate the 
NB. Ohman effect polarization from the Doppler shift. 
NB. Usage Z=. 'atmos_id' spdM lamc tau0 w
NB. A=. ('marcs_package';'p4000_g+3.5') spdME 4000 10 2.5
load'/your_path/Set_LMN.ijs'
NB. Set_LMN'' must be run first: 
Set_LMN''
spdME=: dyad define
'lamc tau0 width'=: y
NB. solve for M-E line in MARCS continuum at frequency points "xx"
'mu0 lams Flxs III QQQ'=: x Extract_line lamc, tau0 
wd 'msgs'                  NB. clear print buffer
th=. cdtr*th=. i.181       NB. 181 theta points, [0,pi]
phi=. (1p1%180)*1.5*i: 60  NB. 121 phi points, -pi/2 < phi < pi/2
'mu xi'=. Xi th; phi
Dshift=: width*(sin phi)*/sin th  NB. Doppler shift of each point
EZ=: th;phi;mu;xi                 NB. make available for inspection
NB. 'th phi mu xi'=. EZ
'I Q'=. Atmos_Pol mu
Q=. Q* cos +:xi  NB. rotate Stokes to z-axis
U=. Q* sin +:xi
dA=. (sin BIN th)*DEL th
dA=. mup*(($mup=.BN mu)$dA)
II=: 1p_1* +/(DEL phi)*(BIN +/"2 dA*BIN"2 I)
QQ=: 1p_1* +/(DEL phi)*(BIN +/"2 dA*BIN"2 Q)
UU=: 1p_1* +/(DEL phi)*(BIN +/"2 dA*BIN"2 U)
POL=: 100*QQ%II
)

NB. Xi is angle between projected surface normal and
NB. Z', the projected Z-axis.
NB. Usage: Xi th;ph --> mu;xi
Xi=: monad define
'th ph'=. y
n=. ph ,&#th              NB. dimensions: phi x del
mu=. (cos ph)*/sin th     NB. ph columns, delta rows    
xi=. atan (sin ph)*/tan th
mu; xi
)

NB. Interpolate I and Q for all values of mu.
Atmos_Pol=: monad define
nn=. #mur=. ,y              NB. ravel mu's to single dimension
IT=: QT=:  (nn,(#lams))$0   NB. empty 2D array to be filled below
DI=. lams splineA III       NB. derivatives for spline interpolation of I
DQ=. lams splineA QQQ       NB. derivatives for spline interpolation of Q
i=. _1
while. (i=. >:i)<nn do.
  shifted_lams=. lams+ i{,Dshift       NB. lambda + Dshift
  IS=: (lams;III;DI) spltrpA shifted_lams
  IS=: (mu0;IS) splint i{mur           NB. interpolate at specific mu
  QS=: (lams;QQQ;DQ) spltrpA shifted_lams      
  QS=: (mu0;QS) splint i{mur
  IT=: IS i} IT                        NB. put shifted profile into array
  QT=: QS i} QT
end.
IT=: (($y),(#lams))$,IT                NB. reformat back to theta,phi array
IT;(($y),(#lams))$,QT
)

NB. B=: Bnu&nu0=. 3e18%x   NB. nu = c%lambda
BN=: monad :'BIN"1 y'      NB. BIN over theta axis
cdtr=: (1p1%180)           NB. for degree->radians
fold0=: ],~ [: -[: |.}.
fold=: ],~ [: |.}.

NB. Polarization across Milne-Eddington absorption line.
NB. Z=. ('marcs_package';'p3000_g+4.0',tag) Extract_line lam,tau0 
NB. where lam is *the* central wavelength of the line and tau0
NB. the ratio of kappa_line/kappa_cont at line center.
NB. Z=. ('marcs_package';'sun.opa') Extract lam
NB. The output is: 'mu rad deg xx Flxs Ia Qa pol'=. Z  
NB. Where the 1st 3 are the angles: mu theta(rad) theta(deg)
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/splines.ijs'
load'/your_path/lintrp.ijs'
load'/your_path/IQMsi.ijs'
load'/your_path/IQs.ijs'
Extract_line=: dyad define
'lamc tau0'=: y
nu=: 3e18% lamc
'fd fn'=: x   NB. fd=file directory, fn=temp&g part pf file name
dir=: '/your_path/',fd,'/'
Z=: mread dir,(fn,tag)  NB. read in the .opa file
AA=. Z do_file lamc    NB. <== extract data from file format           
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  
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. form partial coefficient matrix QQ=[a,b/c,d]
frame=. < (3,(#mua))$mua, theta, (180*theta%1p1)
'tau T ops absc sca'=: AA      NB. pull off data for lamc
T=. T,~ (tau;T) splint 0
ops=. ops,~ (tau;ops) splint 0
absc=. absc,~ (tau;absc) splint 0
sca=. sca,~ (tau;sca) splint 0
kapl=. tau0* ^- *:xx=. 0.1*i. 81    NB. frequencies out to +/- 8 Doppler widths
tau=. 0, tau
ind=. _1
while. (ind=. >:ind)<81 do.         NB. loop over line profile
 abs=. (1+ind{kapl)*absc            NB. absorption = 1+ tau0 x continuum
 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)
 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
 's p'=: QQ s_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'=: (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)+(+/PHI4*|:frame_p)
print 'Wavelength:';lamc;'   Flux:'; {.Flxs    
Z=. >frame
'mu rad deg'=. 0{Z
'Ia Qa pol'=. |:"2|:|:"2 z1=. }.Z
Ia=. |: fold Ia
Qa=. |: fold Qa
mu;(fold0 xx);(fold Flxs);Ia;Qa
)

s_p=: dyad define
'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 
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"
)

do_file=: dyad define  NB. extract data from format of file
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
NB. Note that abs and sca are scaled: abs=kappa/ops, sca=sigma/ops !
NB. In MARCS fortran utility to read .opa files, we find the code:
NB.    abs(i,k)=abs(i,k)*ops(k)
NB.    sca(i,k)=sca(i,k)*ops(k)
NB. which is applied to abs ans scs before they are output.
Qa=. (lam;|:abs) lintrp y 
Qs=. (lam;|:sca) lintrp y 
(A,Qa),Qs             NB. append abs(lam) & sca(lam) to each set
)

mread=: (0&".)@('m'&fread)  NB. read in data matrix
Bnu=: 1.47449934e_47&*@(^&3@])% -&1@(^@(4.799243e_11&*@(] % [)))
NB. last part of name for models in "marcs_package":
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'