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 array
NB. and inc the angle of inclination (degrees).
NB. Usage: Z=. lambda spin_pol_atmos 25700 0.95 90
NB. Then do  'lam I Q U Pol Chi'=. Z
load '/Users/jph/J6/brent.ijs'
spin_pol_atmos=: dyad define
mm=: # lambda=: x
print lambda
'T w inc'=: y
print T,w,inc
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
EZ=: th;phi;delta;n_dot_r;xx;gg;Tt;mu;xi 
NB. 'th phi delta n_dot_r xx gg Tt mu xi'=. EZ
'IA QA'=. Atmos_Pol Tt;mu
IIA=.QQA=.UUA=.POLA=.chiA=. mm$0
k=. _1
while. (k=. >:k) < mm do.  NB. k index loops over lambdas
I=. k{IA
Q=. k{QA
print ' dimension of Q:  ',(": $Q),',  wavelength: ',(": k{lambda)
Q=. Q* cos +:xi  NB. rotate Stokes to z-axis
U=. Q* sin +:xi
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)
UU=. 2p_1* +/(DEL phi)*(BIN +/"1 dA*BN U)
POL=. 100*QQ%II
chi=. -: (180%1p1)*atan2 UU,QQ
print 'Lambda, I, Q, U, Pol, Chi'
print (": (k{lambda),II,QQ,UU,POL,chi)
IIA=. II k}IIA
QQA=. QQ k}QQA
UUA=. UU k}UUA
POLA=. POL k}POLA
chiA=. chi k}chiA
end.
print '------------------------------------------------'
lambda;IIA;QQA;UUA;POLA;chiA
)

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
mask=. (y=0)+.(y=1p1)   NB. fix poles
(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.18)*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
T_p=: Tnew
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
)

NB. Interpolate I and Q for all values of mu.
Atmos_Pol=: monad define
'Tt mu'=. y     NB. Tt: temps for every latitude th (1801)
nn=: $mu        NB. nn <-- (#phi),(#theta) e.g. (901,1801)
mm=: #lambda
FLX=: (mm,(#Tt))$0
T0=. 20000 23000 27500
FLX=: (T0;Flx0) lintrp Tt
print 'dimension of FLX:  ',(": $FLX)
ITT=. QTT=. (mm,nn)$0
k=. _1
while. (k=. >:k) < mm do.  NB. k index loops over lambdas
IxI=: k{"1 III             NB. I(mu0) @ 3 T0 --> 3,159
IoT=: (T0;IxI) lintrp Tt   NB. I(mu0) for every Tt(theta) -> 1801,159
QxQ=. k{"1 QQQ
QoT=. (T0;QxQ) lintrp Tt  NB. dimensions (#theta),(#mu0)
IT=: QT=. nn$0
i=. _1
while. (i=. >:i)<(1801) do.  NB. loop over all Tt latitudes
   ll=. _1
   while. (ll=. >:ll)<(901) do. NB. loop over all longitudes
  a=: (mu0;(i{IoT)) lintrp (i{ll{mu)  NB. mu0 fixed mu array 
  b=. (mu0;(i{QoT)) lintrp (i{ll{mu)
  IT=: a (<ll,i)} IT
  QT=. b (<ll,i)} QT
end. 
end.
ITT=. IT k}ITT
QTT=. QT k}QTT
end.
print 'dimensions of ITT:  ',(": $ITT)
ITT;QTT
)

NB. Bnu=: 1.47449934e_47&*@(^&3@])% -&1@(^@(4.799243e_11&*@(] % [)))
NB. B=: Bnu&nu0=. 3e18%x   NB. nu = c%lambda
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 lambdas.
read_atmos_pol=: monad define
load '/Users/jph/J6/SPLINES/lintrp.ijs'
load '/Users/jph/J6/SPLINES/splines.ijs'
mu0=: ,mread 'atmos_mu'
lams=: ,mread 'atmos_lams'
Flx20368=. ,mread 'Flx_20_3.68'
Flx23393=. ,mread 'Flx_23_3.93'
Flx27424=. ,mread 'Flx_275_4.24'
I20368=. mread 'I_20_3.68'
I23393=. mread 'I_23_3.93'
I27424=. mread 'I_275_4.24'
Q20368=. mread 'Q_20_3.68'
Q23393=. mread 'Q_23_3.93'
Q27424=. mread 'Q_275_4.24'
D=. lams spline Flx20368
Flx20368=. (lams;Flx20368;D) spltrpA lambda=. y
D=. lams spline Flx23393
Flx23393=. (lams;Flx23393;D) spltrpA lambda
D=. lams spline Flx27424
Flx27424=. (lams;Flx27424;D) spltrpA lambda
D=. lams splineA |:I20368
I20368=. (lams;(|:I20368);D) spltrpA 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
Flx0=: >Flx20368;Flx23393;Flx27424
)

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
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'