NB. Polarization across line in a slowly rotating (spherical) 
NB. star using MARCS model atmosphere tables. Lamb the  
NB. wavelength in A. The inclination is fixed at 0 deg.
NB. This uses a Milne-Eddington line profile to evaluate the 
NB. Ohman effect polarization from the Doppler shift. 
NB. Uses radial-tangential model for macroturbulence.
NB. Usage: Z=. 'atmos_id' spdF lamc tau0 rot mact rt_fract
NB. mact = macroturbulent line width; rot = rotational velocity
NB. E.g: Z=. ('marcs_package';'p4000_g+3.5') spdF 4000 10 2 3 0.5
NB. 'line lams Flxs I Q U P'=. Z ;line=lamc,tau0,rot,mact,rt_fract
load'/your_path/Set_LMN.ijs'
NB. Load the "J" calls to the fast Fourier transform libraries:
load'math/fftw'     NB. loads verbs "fftw" and inverse "ifftw"
spdF=: dyad define
Set_LMN''
'lamc tau0 rot mact rt_fract'=: 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=: rot*(sin phi)*/sin th  NB. rot = max rotational Doppler shift           
'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=.BIN"1 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. QQ divided by max value of II
(x;y);lams;Flxs;II;QQ;UU;POL
)

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 rows, th columns       
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 grid to single dimension (y=mu)
DI=. mu0 splineA Iinv=. |:III  NB. derivatives for spline interpolation in mu
DQ=. mu0 splineA Qinv=. |:QQQ     
IT=. QT=. (nn,(#lams))$0       NB. empty 2D array to be filled below
i=. _1
while. (i=. >:i)<nn do.
  shifted_lams=. lams+ i{,Dshift       NB. lambda + Dshift
  mui=. i{mur                          NB. mu value for index "i"
  IS=. (mu0;Iinv;DI) spltrpA mui       NB. profile for local mu
  IB=. IS FTmacro mui                  NB. I-profile broadened by macroturbulence
  IS=. (lams;IB) splint shifted_lams   NB. Doppler shift the profile
  QS=. (mu0;Qinv;DQ) spltrpA mui    
  QB=. QS FTmacro mui                  NB. Q-profile broadened by macroturbulence 
  QS=. (lams;QB) splint shifted_lams 
  IT=. IS i} IT                        NB. put shifted profiles into array
  QT=. QS i} QT
end.
IT=. (($y),(#lams))$,IT                NB. reformat back to theta,phi array
IT;QT=. (($y),(#lams))$,QT
)

NB. IN FTmacro mu -> convolves IN with R/T Gaussians with width mact
FTmacro=: dyad define
INF=. fftw x              NB. Fourier transform of IN (I or Q profile)
ConR=. (x;INF) Rprofile y
ConT=. (x;INF) Tprofile y
NB. multiply radial & tangential components by fractional areas:
(rt_fract*ConR)+(1-rt_fract)*ConT  NB. sum R and T broadened profiles 
)
phiG=: dyad : '(%1p0.5*x)*^-*:y%x'  NB. normalized Gaussian profile
NB. (IN;INF)Rprofile mu -> IN profile broadened by radial motions
Rprofile=: dyad define 
'IN INF'=. x
radial=. mact*y         NB. y = cos theta
if. radial< 0.025 do.   NB. if broadening too small, return IN profile
IN  else.
scale=. 1{DEL lams      NB. size of sample interval
GPrF=. fftw (radial phiG lams)
ConFr=. GPrF * INF
ConR=. Re ifftw ConFr
nn=. i. hn=. >. -: #lams
scale* fold nn{ConR  end.
)
NB. (IN;INF)Tprofile mu -> IN profile broadened by tangential motions
Tprofile=: dyad define 
'IN INF'=. x
tangent=. mact*(%:1-y*y) NB. (%:1-y*y) = sin theta
if. tangent< 0.025 do.   NB. if broadening too small, return IN profile
IN  else.
scale=. 1{DEL lams       NB. size of sample interval
GPtF=. fftw (tangent phiG lams)
ConFt=. GPtF * INF
ConT=. Re ifftw ConFt
nn=. i. hn=. >. -: #lams
scale* fold nn{ConT  end.
)

NB. B=: Bnu&nu0=. 3e18%x   NB. nu = c%lambda
cdtr=: (1p1%180)           NB. for degree->radians
fold0=: ],~ [: -[: |.}.
fold=: ],~ [: |.}.

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 is
NB. the ratio of kappa_line/kappa_cont at line center.
NB. Example: 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) 
NB. Try: 'pensize 3' plot mu;pol
require 'files'
load'/your_path/SPLINES/splines.ijs'
load'/your_path/SPLINES/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=: '/Users/jph/J6/Atmos_Pol.d/',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
nl=. 121                            NB. 12 Doppler widths
kapl=. tau0* ^- *:xx=. 0.1*i. nl
tau=. 0, tau
ind=. _1
while. (ind=. >:ind)<nl do.         NB. loop over line profile
 abs=. (1+ind{kapl)*absc            NB. absorption increased to 10x 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'
max=: >./