NB. calculate polarization of a rapidly rotating star
NB. using model atmosphere tables. T = effective temp
NB. of the non-rotating counterpart, w the fraction
NB. of break-up velocity, lambda the wavelength
NB. and inc the angle of inclination (degrees).
spin_pol_atmos=: monad define
'T lambda w inc'=: y
inc=. cdtr*inc
read_atmos_pol lambda
NB. M=. 4.0   NB. spectral type B5 V (T=15,200 K)
NB. L=. 240   NB. mass and luminosity (solar units)
M=. 10        NB. spectral type B1 V (T=25,700 K)
L=. 5754
g_pole=: 10^. 2.46e_11*(M%L)*T^4
th=. cdtr*th=. 0.1*i.1801  NB. 1801 theta points, [0,pi]
'delta n_dot_r'=. w spin_delta th
'th xx gg Tt'=. T_of_theta T,w
phi=. (1p1%1800)*i. 901
'mu xi'=. inc Xi th; phi; delta
'I Q'=. Atmos_Pol Tt;mu
Q=. Q* cos +:xi  NB. rotate Stokes to z-axis
dA=. (BIN n_dot_r)%~(BIN *:xx)*(sin BIN th)*DEL th
mup=. BN mu
dA=. mup*(($mup)$dA)
II=. 2p_1* +/(DEL phi)*(BIN +/"1 dA*BN I)
QQ=. 2p_1* +/(DEL phi)*(BIN +/"1 dA*BN Q)
POL=. 100*QQ%II
II,QQ,POL
)

NB. Shape of rotating star x(theta,w) where
NB. w=vel/v_break-up & theta is angle
NB. wrt z-axis in radians.    
NB. Usage: w star_x theta --> radius(theta)
star_x=: dyad define
a=. 3% (w=.x)* sth=. sin y
b=. 3%~ 1p1+ acos w*sth
a* cos b
)

NB. effective gravity of rotating star in units
NB. of GM/R_p^2. Usage: w spin_g theta --> g_eff
spin_g=: dyad define
X=. (w=.x) star_x y
a=. 1- (8%27)*w*w*X^3
b=. %:(*:cos y)+ (*: a*sin y)
b% *:X
)

NB. Delta is the angle delta between the surface normal
NB. and the z-axis (rotation axis) of the star. n_dot_r
NB. "n_dot_r" is the cos of the angle between surface normal
NB. and the radial direction from the center.
NB. Usage: w star_delta theta --> delta (radians), n dot r.
spin_delta=: dyad define
X=. (w=.x) star_x (th=. y)
sf=. 1- (8%27)*w*w*X^3
delta=. 1p1|1p1+ atan sf*tan th
NB. Note: delta<th for th<pi/2 & delta>th for th>pi/2
NB. returns delta in 2nd quadrant for theta in 2nd quad.
delta=. 1p1 (I. th=1p1)} delta   NB. fix tan=0 end-point
n_dot_r=. cos th-delta
delta; n_dot_r
)

T_of_theta=: monad define
'T w'=. y
dth=. 0.0001
thr=. cdtr* th=. dth*i. 900001
xx=. w star_x thr
'delta n_dot_r'=. w spin_delta thr
dA=. xx*xx%n_dot_r
dA=. 2p1*dth*dA*sin thr
A=: 360%~ +/dA    NB. ratio of surface to 4 pi.
gg=. w spin_g thr
dL=. gg*dA                    NB. flux per d(theta)
L=: 360%~ +/dL                NB. luminosity vs. 4piR^2 sig T^4
Tnew=: (%L^0.25)*T            NB. new polar (gg=1) temp
dth=. 0.1
thr=. cdtr*th=. dth*i. 1801   NB. use 1801 theta points, 0 to pi
thx=. cdtr*thh=. dth*>:i. 900 NB. use 900 theta points, 0.1 to pi/2
xx=. 1, (w star_x thx)        NB. prepend the thx=0 point
xx=. xx, |. }:xx              NB. pi/2 - pi is reverse of 0 - pi/2
gg=. 1, (w spin_g thx)
gg=. gg, |. }:gg
Tt=. Tnew*(gg^0.25)
>thr;xx;gg;Tt
)

NB. Xi is angle between projected surface normal and
NB. Z', the projected Z-axis.
NB. Usage: inc Xi th;ph;del --> mu;xi
Xi=: dyad define
sinc=. sin inc=. x
cinc=. cos inc
'th ph del'=. y
NB. ((sin a)*/(sin b))+(a ,&# b)$d
n=. ph ,&# del  NB. dimensions: phi x del
mu=. |((cos ph)*/(sin del)*sinc)+ n$ ,cinc*cos del
a=. (sin ph)*/(sin del)    NB. ph columns, delta rows
b=. (n$ ,sinc*cos del)- (cinc*cos ph)*/sin del
NB. need atan2 so xi>p/2 for lower hemisphere:
xi=. |:atan2 |:>a;b
mu; xi
)

NB. Interpolate I and Q for all values of mu.
Atmos_Pol=: monad define
'Tt mu'=. y
T0=. 20000 23000 27500
n=. $mu
IoT=. (T0;III) lintrp Tt  NB. I(mu0) for every Tt(theta)
QoT=. (T0;QQQ) lintrp Tt  NB. dimensions (#theta),(#mu0)
IT=. QT=. |: n$0
mut=. |: mu               NB. transpose of mu array
i=. _1
while. (i=. >:i)<#Tt do.
  a=. (mu0;(i{IoT)) lintrp (i{mut)
  b=. (mu0;(i{QoT)) lintrp (i{mut)
  IT=. a i} IT
  QT=. b i} QT
end.
(|:IT);(|:QT)
)

BN=: monad :'BIN"1 y'      NB. BIN over theta axis
cdtr=: (1p1%180)

NB. Read in atmosphere tables for I,Q and Flux for model
NB. Teff & log g, and interpolate at the requested lambda.
read_atmos_pol=: monad define
load '/your_path/splines.ijs'
mu0=: ,mread '/your_path/atmos_mu'     NB. read mu's for which I&Q tabulated
lams=: ,mread '/your_path/atmos_lams'  NB. read lambdas for which I&Q tabulated
I20368=. mread 'I_20_3.68'             NB. read I(mu) for T=20,000K, log g=3.68 
I23393=. mread 'I_23_3.93'
I27424=. mread 'I_275_4.24'
Q20368=. mread 'Q_20_3.68'             NB. read Q(mu) for T=20,000K, log g=3.68
Q23393=. mread 'Q_23_3.93'
Q27424=. mread 'Q_275_4.24'
D=. lams splineA |:I20368
I20368=. (lams;(|:I20368);D) spltrpA lambda  NB. interpolate I&Q at desired lambda
D=. lams splineA |:I23393
I23393=. (lams;(|:I23393);D) spltrpA lambda
D=. lams splineA |:I27424
I27424=. (lams;(|:I27424);D) spltrpA lambda
D=. lams splineA |:Q20368
Q20368=. (lams;(|:Q20368);D) spltrpA lambda
D=. lams splineA |:Q23393
Q23393=. (lams;(|:Q23393);D) spltrpA lambda
D=. lams splineA |:Q27424
Q27424=. (lams;(|:Q27424);D) spltrpA lambda
III=: >I20368;I23393;I27424
QQQ=: >Q20368;Q23393;Q27424
)