NB. Calculate Polarization from a binary star distorted NB. to a Roche potential surface. Assumes emergent I & Q NB. from a pure electron scattering atmosphere. NB. Usage: Binary q, scale, phase, inclin NB. phase= phase angle; inclin= inclination of z-axis to observer NB. (phase and inclin are input in degrees) NB. Example (input): IN45=. |:>(19#1);(19#1.01);(10*i. 19);(19#45) NB. q=1;Omega=1.01*Omega(critical);phase=1,5,10,...175,180;inc=45 NB. Example (run): ES45=. Binary"1 IN45 NB. Example (output): 'phs inc pA I Q U Pol chi'=. |:ES45 NB. 'pensize 3' plot phs; Pol <-- plot Pol vs phase load'/your_path/splines.ijs' load'brent.ijs' NB. solve equations by Brent's method NB. "brent" is found near the bottom of this webpage cdtr=: (1p1%180) NB. convert degrees to radians Binary=: monad define read_pure_scat'' 'q scale phase inclin'=: y xcm=: q% q1=: 1+ q NB. xcm= center of mass x-coordinate qi=: 2%q1 phi0=: cdtr*phase NB. phi0 = phi of viewer's phase inc=: cdtr*inclin rL1=: RL1 q OM1=. OML1 q OM=: OM1 * scale NB. scale > or = 1, choose potential do_thph'' ) do_thph=: monad define theta=: 10* i. nth=. 19 NB. 19 theta points, [0,180] nph=. # phi=: 10* i: 18 NB. 37 phi points, [-180,180] theta=: 5* i. nth=. 37 NB. 37 theta points, [0,180] nph=. # phi=: 5* i: 36 NB. 73 phi points, [-180,180] theta=: i. nth=. 181 NB. 181 theta points, [0,180] nph=. # phi=: i: 180 NB. 361 phi points, [-180,180] NB. theta=. -: i. nth=. 361 NB. 361 theta points, [0,pi] NB. nph=. # phi=. -: i: 360 NB. 721 phi points, [-pi,pi] nn=: */ dim=: nth,nph thA=: |:(|.dim)$cdtr*theta phA=: dim$cdtr*phi gg=:RR=: nn#0 NB. create empty arrays Ndc=:Grav=: (nn,3)$0 ii=. _1 NB. ---------------------------------------------------- while. (ii=. >:ii)0 do. rmin=: OM_rgrad brent (0.1,0.8,1e_10,0) end. r=. OM_eval brent (0.1,rmin,2e_10,OM) RR=: r ii}RR nx=. (cos ph)*sin th ny=. (sin ph)*sin th nz=. cos th dc=. nx,ny,nz NB. direction cosigns to th,ph Ndc=: dc ii}Ndc NB. triplet of direction cosigns 'rx ry rz'=. r*dc NB. coordinates of surface point r13=. 3^~ r1=. %:(rx^2)+(ry^2)+(rz^2) NB. distance to M1 r23=. 3^~ r2=. %:((1-rx)^2)+(ry^2)+(rz^2) NB. distance to M2 a=. (-qi*rx%r13)+(2*xcm*(1-rx)%r23) OMx=. a+ 2*(rx-xcm) a=. (-qi*ry%r13)+(-2*xcm*ry)%r23 OMy=. a+ 2*ry OMz=. (-qi*rz%r13)+(-2*xcm*rz)%r23 Grav=: (OMx,OMy,OMz) ii}Grav end. NB. --------------------------------------------------- RR=: dim$RR NB. r of each surface point Rxyz=: RR* Ndc=: dim$Ndc NB. grav=: L Grav=: dim$Grav Norm=: -Grav% grav NB. -g/|g|= surface normal gamma=: +/"1 Ndc*Norm NB. cos of angle between r & -g RNorm=. Norm qRot3 ((0 0 1);-phi0) NB. rotate about z-axis by phi0 RNorm=. RNorm qRot3((0 1 0);(0.5p1-inc)) NB. about y'-axis by pi/2-inc mu=: 0{"1 RNorm NB. x component = mu viewed from (phi0,inc) dA=. RR*RR%gamma NB. r^2/(n dot r) dA=. 2p1*dA*sin thA A=: (180*360)%~ +/+/dA NB. total surface area RNy=: 1{"1 RNorm RNz=: 2{"1 RNorm xi=: atan2 |:>(,RNy);(,RNz) NB. angle between n & z-axis sign=: _2p1*(xi>:1p1) NB. sign= -2pi if xi >= pi xi=: dim$ sign+xi NB. map xi>=p1 into xi<0 values 'Ia Qa'=: grav*"2 >Interp mu NB. Set Ia & Qa proportional to g Q=: Qa* cos +:xi NB. rotate Stokes to z-axis U=: Qa* sin +:xi Dph=: {. DEL"1 phA mup=. BIN mu*(mu>0) dAA=: (BIN gamma)%~mup*(BIN *:RR)*(sin BIN thA)*(DEL thA) pA=. 2p_1* +/Dph*BIN(+/dAA) NB. pA = projected area of star II=. 2p_1* +/Dph*BIN(+/dAA*BIN Ia) QQ=. 2p_1* +/Dph*BIN(+/dAA*BIN Q) UU=. 2p_1* +/Dph*BIN(+/dAA*BIN U) POL=. 100* (%:(*:QQ)+(*:UU))%II chi=. -: (180%1p1)*atan2 UU,QQ phase,inclin,pA,II,QQ,UU,POL,chi ) NB. 'surface' plot theta;phi;mu NB. Calculate Roche potential as function of r NB. th, ph & q set externally as global nouns NB. Usage: OM_eval r --->OM OM_eval=: monad define xcm=. q% q1=. 1+ q sth=. sin th cph=. cos ph trr=. y* lam=. sth*cph b=. % %: 1+(*:y)- 2*trr ((2%q1)%y)+((2*xcm)*b-trr)+(*:y*sth)+ *:xcm ) NB. Calculate Roche potential as function of r NB. th, ph & q set externally as global nouns NB. Usage: OM_rgrad r --->gradient OM_rgrad=: monad define r=. y sth=. sin th cph=. cos ph xcm=. q% q1=. 1+ q qi=. 2%q1 trr=. r* lam=. sth*cph c=. a* b=. %: a=. % 1+(*:r)- 2*trr (-(2%q1)%r*r)+(2*xcm*((lam-r)*c)-lam)+ 2*r* *:sth ) NB. get x coordinate of 1st Lagrange point NB. for given q=M2/M1. Usage RL1 q ---> r(L1) NB. Solve equation dOM/dx = 0 for case y=z=0. RL1=: monad define y&fss brent (0.05,0.95,1e_12,0) ) fss=: 4 : '(-%*:y)+(x% *:(1-y))+((1+x)*y)-x' NB. get potential of L1 surface for given q OML1=: monad define L1=. RL1 q=. y xcm=. q% q1=. 1+ q ((2%q1)%L1)+(2*xcm%(1-L1))+(L1-xcm)^2 ) NB. Rotate n vectors v about axes u by angles w. NB. Makes use of this quaternion-based formula: NB. with u'= sin(w/2)*u and t= 2(u' X v) NB. v' = v + [t * cos(w/2)] + (u' X t) NB. Dimensions of arrays are v and u $V --> (n 3). NB. The axes u must be normalized to unit length. NB. Usage: v qRot3 (u;w) --> v' (v' rotated v's) NB. Any of v,u & w may be single values or arrays. qRot3=: dyad define 'u w'=. y NB. u=. norm u NB. Add this if u's are not normalized. if. 2>#$u do. u=. ($x)$u end. NB. for single axis, clone it if. (2>#$u) *. (1<#w) do. NB. for single u & v ... u=. u */~ sin-:w NB. this will produce (#w) u's else. u=. u * sin-:w NB. if there's a u for each w, "*" OK end. NB. for one axis and one vector v, but multiple w's, clone v: if. 2>#$x do. x=. ($u)$x end. t=. +: u cross x x+ (t*cos-:w)+ u cross t ) L=: +/"1 &. (*:"_) NB. Lengths of vectors norm=: ] % L NB. normalize to unit length NB. cross products of vectors: v1 cross v2 --> v1 X v2 cross=: (([: 1&|."1 [)* [: 2&|."1 ])- (([: 1&|."1 [)* [: 2&|."1 ])~ L=: +/"1 &.(*:"_) fold0=: ],~ [: -[: |.}. fold=: ],~ [: |.}. NB. Interpolate I and Q for all values of mu. Ravel 2D mu array NB. and restore shape after interpolation. Interp=: monad define mu=. ,y mask=. mu>0 NB. mask=0 where mu <= 0 I=. (mu0;I_es;0{Dspline) spltrp mu Q=. (mu0;Q_es;1{Dspline) spltrp mu (dim$mask*I);(dim$mask*Q) ) NB. Use the pure scattering solution (Chandrasekhar, "Radadiative NB. Transfer", p248). 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 Q_es=: -pol*I_es NB. Q<0 as pol perp to normal Dspline=: mu0 splineA >I_es;Q_es ) NB. by J. Patrick Harrington, 27 Nov 2016