NB. To calculate polarization of a rapidly rotating star
NB. with a pure electron scattering atmosphere. T is the
NB. effective temperature of the non-rotating counterpart,
NB. w is the fraction of break-up velocity, lambda the 
NB. wavelength and inc the angle of inclination (degrees).
NB. This version uses Espinosa Lara & Rieutord (A&A 533, A43)
NB. gravity darkening theory rather than von Zeipel's law.
load '/your_path/brent.ijs'
spin_pol_es2=: monad define
read_pure_scat''
'T w inc'=: y
inc=. (1p1%180)*inc
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;gg;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=omega/omega_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
mask=. (y=0)+.(y=1p1) NB. fix pole singularities
(mask* (#y)#1)+(-. mask)* 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)
xx=. xx, |. }:xx
gg=. 1, (w spin_g thx)
gg=. gg, |. }:gg
T_p=: Tnew
NB. Solve Esponisa Lara & Rieutord equation for T(theta)
Tt=. ,(w&Tvar"(0)(0,thx))
Tt=. Tt, |. }: Tt
>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
)

Atmos_Pol=: monad define
'Tt gg mu'=. y
Fmu=. (Tt%T)^4      NB. flux varies with "th" axis: Tt(theta)
'I Q'=. Interp mu
(Fmu*"1 I);(Fmu*"1 Q)  NB. scale by sig*T^4
)

BN=: monad :'BIN"1 y'
cdtr=: 1p1%180         NB. convert degrees to radians

NB. Interpolate I and Q for all values of mu. Ravel 2D mu array
NB. and restore shape after interpolation.
Interp=: monad define
n=. $mu=. y
I=. n$ (mu0;I_es;0{Dspline) spltrp ,mu
Q=. n$ (mu0;Q_es;1{Dspline) spltrp ,mu
I;Q
)

NB. Use the pure scattering solution (Chandrasekhar, "Radadiative
NB. Transfer", p 248). Compute derivatives for spline interpolation.
NB. Check: +: mu0 sinteg I_es*mu0 --> 1.00002 
read_pure_scat=: monad define
mu0=: 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 
mu0=: mu0, 0.7 0.75 0.8 0.85 0.9 0.95 1
I_es=: 0.41441 0.47482 0.52397 0.57001 0.61439 0.65771 0.70029 0.74234 
I_es=: I_es, 0.78398 0.8253 0.86637 0.90722 0.94789 0.98842 1.02882 
I_es=: I_es, 1.06911 1.10931 1.14943 1.18947 1.22945 1.26938
pol=. 0.11713 0.08985 0.07448 0.06311 0.0541 0.04667 0.04041 0.03502 
pol=. pol, 0.03033 0.02619 0.02252 0.01923 0.01627 0.01358 0.011124 
pol=. pol, 0.00888 0.006822 0.004919 0.003158 0.001522 0
load '/your_path/SPLINES/splines.ijs'
Q_es=: -pol*I_es   NB. Q<0 as pol perp to normal
Dspline=: mu0 splineA >I_es;Q_es
)

NB. w7=. 0.925723   NB. w for om=0.7
NB. w9=. 0.9928072  NB. w for om=0.9
NB. w5=. 0.76980036 NB. w for om=0.5
L=. 335      NB. Parameters for Regulus
R_pole=. 3.22
T_p=: 5779.19*(L % R_pole^2)^(%4)

NB. Temperature law from 2011 paper by
NB. Espinosa Lara & Rieutord, A&A 533, A43.
NB. Usage: w Tvar theta(radians) --> T(th)
Tvar=: 4 : 0
om=: w2om w=. x   NB. find omega corresponding to w
th=. y
xeq=. w star_x 0.5p1
rt=. xeq%~ w star_x th  NB. rt=r_tilde
  a=. (3%~ om*om)*(rt^3)*(cos th)^3
  b=. (cos th)+ ^. | tan -:th
  tau=: a+b             NB. tau as in eq (24)
if. y=0 do. 
  A=. 0.25^~ ^(2%3)*om*om*rt^3    NB. eq (27)
elseif. y=0.5p1 do. 
  A=. 0.25^~ (1-om*om)^(_2%3)     NB. eq (28)
elseif. y>0 *. y<0.5p1 do.
  TH=. ff brent 0,1.6,1e_8,tau    NB. eq (24)
  A=. %: (tan TH)%(tan th)
end.
b=. (%rt^4)+ (om^4)*rt*rt*sin th
b=. b- (2%rt)*(om* sin th)^2
T=. T_p*(A% %:xeq)* b^(%8)       NB. eq (31)
)

ff=: 3 : '(cos y)+^. tan -:y'
NB. w2om w --> om = Omeg/Omeg_Kepler
w2om=: 3 : 'y*(1.5%~ y star_x 0.5p1)^1.5'