NB. get light curve & polarization of transits
NB. Usage: out=. (mua;Ia;Qa) transit (Rs, Rp, ip)
NB. where (mua;Ia;Qa) are the stellar atmosphere parameters:  
NB. mua = array of mu=cos(theta)'s, Ia=I(mu), Qa=Stokes Q(mu) 
NB. Rs= star radius, Rp= planet radius, ip= impact parameter
NB. (closest approach to star center) of transit track
NB. 1st row of "out" is 201 equally spaced points on transit 
NB. track, 2nd row = distance from star center, 3rd = eclipse 
NB. depth, 4th row = fractional polarization of starlight                  
transit=: dyad define
'mua Ia Qa'=. x
'Rs Rp ip'=: y
NB. Total from star outside transit:
Itot=: 1p1*(*:Rs)*(+: +/(BIN Ia*mua)*DEL mua)
ip2=: *: ip
zmax=: %: (*:Rs+Rp)-ip2
z=. zmax* int01 <: nz=. 101
s0=. %: ip2+ *:z
rhs=. |:(2,nz)$ s0,z
fold |: x&disk"(1) rhs
)

load'/your_path/splines.ijs'

NB. Do integrals over s, where s = projected radius from
NB. disk center and s0 = center of transiting object 
disk=: dyad define
'mua Ia Qa'=. x
's0 z'=. y
smax=. Rs <. s0+ Rp
dels=. smax - smin=. |s0- Rp
s=. smin+ dels*int01 1000
'mu alpha arc warc'=. geometry s0;s
I=. (mua;Ia) splint mu
Q=. (mua;Qa) splint mu
Ist=. +/(BIN arc*I)*DEL s
Qst=. +/(BIN warc*Q)*DEL s
if. 0< rc=. (Rp-s0) do.
 s=. rc*int01 300
 mu=. %: 1- *:s%Rs
 I=. (mua;Ia) splint mu
 Icc=. 2p1* +/(BIN s*I)*DEL s
else.
 Icc=. 0
end.
pol=. Qst%(Itot-(Ist+Icc))
depth=. 1- (Ist+Icc)%Itot
z,s0,depth,pol
)
  
NB. geometry of transiting object 
geometry=: monad define
's0 s'=. y
mu=. %: 1- *:s%Rs
s02=. *:s0
r2=. *:Rp
cosa=. (2*s*s0)%~ (s02-r2)+ *:s
alpha=. acos(*cosa)* 1<.|cosa
arc=. 2*s*alpha
warc=. arc* sinc +:alpha
mu;alpha;arc;warc
)

sinc=: 3 : '(y=0)+y%~ sin y'  NB. sinc(x) function

NB. fold data about 1st point, while negating 1st half of 
NB. 0th (z) row, i.e., [0 1 2 3]->[-3 -2 -1 0 1 2 3]:
NB. fold=: monad define '((_1(0)}(#y)#1)*|."1 y),"1(}."1 y)'
fold=: (((_1)0}1 #~ #)*|."1),"1 }."1