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