NB. Functions for solution by matrix approximation to integral transforms.
NB. This code is in the "J" language.

NB. Direct solution for s(tau) & p(tau)
NB. where B(tau) and lambda(tau) are supplied
s_and_p=: dyad define
'L0 M0 N0'=. x         NB. the transform matrices
'lam B'=. y            NB. Plank term B(tau) and lambda(tau)
n=. #L0                NB. n is the number of depth points
idn=. =@i. n           NB. idn is n x n identity matrix
a=. idn- L0* oml=. 1-lam
b=. _3%~oml*M0
c=. (_3%8)*oml*M0
d=. idn - (3%8)*oml*N0
Q=. (a,"1 b),"2 c,"1 d   NB. form grand coefficient matrix  Q=[a,b/c,d]
rhs=. (lam*(n$B)),(n#0)  NB. rhs is the right hand side of the system
sp=. rhs %. Q            NB. solve linear system: rhs %. (coeff matrix)
's p'=. (2,n)$ sp        NB. 1st n elements are "s"; 2nd n are "p"
)
NB. ------------------------------------------------------------------
tau_grid=: monad define
NB. set up log tau grid over 0.001 to 10: 16 points/decade:
ttt=. 2}. log_tau_grid 0.001 16 10
NB. equal tau steps of 1.5 from 10 to 22:
ttt=. ttt,11.5,13,14.5,16,17.5,19,20.5,22
NB. equal tau steps of 0.000143 from 0 to 0.001:
ttt=. (0.001*int01 7),ttt
)
NB. tau1   NB. first optical depth point
NB. npd    NB. number of points/decade
NB. T      NB. Half-depth of slab
NB. SIDE EFFECTS! asigns global variable "nt"
log_tau_grid=: monad define
'tau1 npd T'=. y
t=. T,~(i. +/t<T){ t=. 0,tau1*10^(npd %~ i. 500)  
)
NB. ------------------------------------------------------------------------
Set_LMN=: monad define
NB. Set up tau grid of 80 points for semi-infinite atmosphere:
ttt=: tau_grid ''
NB. evaluate the 80x80 matrix transforms for this tau grid:
L0=: 1 MTsi ttt
M0=: 3 MTsi ttt
N0=: 5 MTsi ttt
PHI2=: {. 2 MTsi ttt
PHI4=: {. 4 MTsi ttt
'      done!  '
)
NB.--------------------------------------------------------------------------
NB. Matrix Transforms for semi-infinite atmosphere.
NB. linear extrapolation from last point to infinity.
NB. idx=1,2,3,4,5 -> Lambda,Phi,M,Phi^4,N
NB. idx MT tau_grid -> (linear);(spline) matrix operator
MTsi=: dyad define
pm=. }:"(1) _1+ +: t<:/t=. y
'I0 I1 I2 I3 IN0 INs'=. x IIsi t
t3=. t* t2=. t*t
a2=. *: a=. }: t
b2=. *: b=. }. t
h2=. *: h=. b-a
h6=. 6%~ hi=. % h
X=. hi*"1 ((-t-/b)*I0)-pm*I1
Y=. -hi*"1 ((-t-/a)*I0)-pm*I1
lin=. (X,"(1) 0)+(0,"(1) Y)
e2=. -({:hi)*INs
e1=. IN0+ ({:hi)*INs
eN=. |: e1(N-1)}e2(N-2)}(N,(N=.#t))$0
lin=. lin+ eN       NB. linear approximation
NB. unit vectors of length rows, columns:
uc=. }. ur=. N#1 
u1=. ur*/ b*(b2-h2)
u1=. u1- t*/((3*b2)-h2)
u1=. u1+ ((3*t2)*/b)-(t3*/uc)
U=. I0*u1
u1=. -ur*/ (3*b2)-h2
u1=. u1+ (6*t*/b)-(3*t2)*/uc
U=. U+ pm*I1*u1
U=. U+ I2*3*(-t-/b)
U=. h6*"1 U- pm*I3
v1=. ur*/ a*(a2-h2)
v1=. v1- t*/((3*a2)-h2)
v1=. v1+ (3*t2*/a)-(t3*/uc)
V=. -I0*v1 
v1=. ur*/ (3*a2)-h2
v1=. v1+ ((3*t2)*/uc)-6*t*/a
V=. V+ pm*I1*v1
V=. V+ I2*3*(t-/a)
V=. h6*"1 V+ pm*I3
C=. SpMD t
spl=. (U mm }:C)+(V mm }.C)
e3=. 6%~({:h)* INs
spl=. spl+ e3*/(N-2){C  NB. spline correction
lin+spl
)
NB. ------------------------------------------------------------------------
NB. make table of I-integrals for MTsi           
IIsi=: dyad define
idx=. x
a=. | y -/ (}:y) NB. |tau(k)-tau(i)|
b=. | y -/ (}.y) NB. |tau(k+1)-tau(i)|
pm=. }:"(1) _1+ +: y<:/y
if. -.(idx=2)+:(idx=4) do. pm=. |pm end.
if. idx=1 do.
NB. the kernels for Lambda transform
K0=. 3 : '_2%~ 2&EIn y'
K1=. 3 : '2%~ (3&EIn y)-(^-y)'
K2=. 3 : '-(4&EIn y)+ 2%~y*(^-y)'
K3=. 3 : '(3*5&EIn y)- 2%~(^-y)*(3+y+y*y)'
end.
if. idx=2 do.
NB. the kernels for Phi-transform
K0=. 3 : '_2*3&EIn y'
K1=. 3 : '+:(2*4&EIn y)- ^-y'
K2=. 3 : '(_12*5&EIn y)+2*(^-y)*(1-y)'
K3=. 3 : '(48*6&EIn y)-2*(^-y)*(6+y*y)'
end.
if. idx=3 do.
NB. the kernels for M-transform
K0=. 3 : '(_2%~ 2&EIn y)+(3%2)*4&EIn y'
K1=. 3 : '-:(3&EIn y)+(_9*5&EIn y)+2*(^-y)'
K2=. 3 : '(-(4&EIn y))+(18*6&EIn y)-(^-y)*(3-y)'
K3=. 3 : '(3*5&EIn y)+(_90*7&EIn y)+(^-y)*(15+(_2*y)+y*y)'
end.
if. idx=4 do.
NB. the kernels for Phi^(4)
K0=. 3 : '((_2%3)*3&EIn y)+(2*5&EIn y)'
K1=. 3 : '((4%3)*4&EIn y)+(_8*6&EIn y)+(4%3)*(^-y)'
K2=. 3 : '(_4*5&EIn y)+(40*7&EIn y)-(4%3)*(^-y)*(4-y)'
K3=. 3 : '(16*6&EIn y)+(_240*8&EIn y)+(4%3)*(^-y)*(24+(_3*y)+y*y)'
end.
if. idx=5 do.
NB. the kernels for N_tau:
K0=. 3 : '((_5%3)*2&EIn y)+(4*4&EIn y)-(3*6&EIn y)'
K1=. 3 : '((5%3)*3&EIn y)+(_12*5&EIn y)+(15*7&EIn y)-(2%3)*(^-y)'
K2=. 3 : '((_10%3)*4&EIn y)+(48*6&EIn y)+(_90*8&EIn y)+(2%3)*(^-y)*(6-y)'
K3=. 3 : '(10*5&EIn y)+(_240*7&EIn y)+(630*9&EIn y)-(2%3)*(^-y)*(63+(_5*y)+y*y)'
end.
I0=. pm*(K0 b) -K0 a
I1=. pm*(K1 b) -K1 a
I2=. pm*(K2 b) -K2 a
I3=. pm*(K3 b) -K3 a
NB. integrals from tau(N) to infinity:
IN0=. -K0 tiN=. ({:y)-y  NB. tiN = tau(N)-tau(i)
INs=. (-K1 tiN)- tiN*IN0 
I0;I1;I2;I3;IN0;INs
)
NB. ------------------------------------------------------------------------
NB. Spline 2nd Derivative Matrix (natural splines): 
NB. (SpMD x) on vector x yeilds matrix B where (B mm y)
NB. gives the spline 2nd derivatives, i.e.,  D = x spline y 
NB. 
SpMD=: monad define
k1=. 6% h1=. }: h=. DEL y
k3=. 6% h3=. }. h
h2=. +: h3+h1
h1=. 0, (}. h1)
h3=. (}: h3),0
AI=. trivert h1;h2;h3      NB. invert (N-2)x(N-2) matrix
k1=. 0,k1,0
k3=. 0,k3,0
k2=. -(k1 + k3)
NB. Expand the three diagonal bands of a  
NB. tridiagonal matrix into the full matrix. 
NB. Rows are  k1, k2, k3
B=. k2* I=. =@i.(#k2)
B=. B+ 1 |."1 (k1* I)
B=. B+ _1 |."1 (k3* I) NB. full N x N matrix
B=. AI mm }: }. B  NB. result is (N-2) x N matrix
0,B,0  NB. result is N x N matrix
)
NB. Invert a tridiagonal matrix by direct substitution.
NB. Rows are  a_i*x_(i-1), b_i*x_i, c_i*x_(i+1)  
NB. a,b,c have same length, i.e., a[0]=0 and c[n]=0. 
NB. Usage: inverse(full nxn matrix)<-- tridgl a;b;c 
NB. 
trivert=: monad define
'a b c'=. y 
n1=. <: n=. <: nn=. #b
x=. %{.b
inv=. x (<0 0)}(=@i. nn)  
v=. (x*{.c) 0}v=. nn#0
j=. 0
while. n1>:j=.>:j do.
  b1=. j{b
  a1=. j{a
  v0=. (<:j){v
  x=. %b1-a1*v0
  c1=. j{c
  v=. (x*c1)j}v
  inv1=. j{inv
  inv0=. (<:j){inv
  inv=. (x*inv1-a1*inv0)j}inv
end.
x=. %(n{b)-(n{a)*(n1{v)
inv=. (x*(n{inv)-(n{a)*(n1{inv)) n}inv
j=. n
while. 0<:j=.<:j do.
  inv1=. j{inv
  v1=. j{v
  inv2=. (>:j){inv
  inv=. (inv1-v1*inv2)j}inv
end.
)
NB. --------------------------------------------------------------------------
NB. Create matrix E to produce emergent intensities from source functions
NB. s(tau) and p(tau) for semi-infinite atmosphere.
IQMsi=: monad define
't mu'=. y
'E0 E1 E2 E3 EN0 ENs'=. mu EEsi t
t3=. t* t2=. t*t
a2=. *: a=. }: t
b2=. *: b=. }. t
h2=. *: h=. b-a
h6=. 6%~ hi=. % h
X=. hi*"1 (b*"1 E0)-E1
X=. X - ({:hi)*ENs
Y=. -hi*"1 (a*"1 E0)-E1
Y=. Y + EN0+ ({:hi)*ENs
lin=. (X,"(1) 0)+(0,"(1) Y)
U=. (b*b2-h2)*"1 E0
U=. U- ((3*b2)-h2)*"1 E1
U=. h6*"1 U+ (3*b*"1 E2)-E3
V=. -(a*a2-h2)*"1 E0 
V=. V+ ((3*a2)-h2)*"1 E1
V=. h6*"1 V- (3*a*"1 E2)-E3
C=. SpMD t
spl=. (U mm }:C)+(V mm }.C)
end=. 6%~({:h)*ENs
spl=. spl+ end*/ {:}:C
lin+spl
)
NB. The 6 E-integrals (= mu^n \int x^n exp(-x) dx) needed by IQMsi 
NB. (I&Q matrix, semi-infinite) to get the emergent rays at angles mu      
NB. Usage: mu EEsi tau --> E0;...
EEsi=: dyad define
a=. }:"1 x%~/y  NB. table of tau_i/mu_j 
b=. }."1 x%~/y  NB. table of tau_(i+1)/mu_j
K0=. 3 : '-(^-y)'
K1=. 3 : '-(^-y)*(1+y)'
K2=. 3 : '-(^-y)*(2+(2*y)+y*y)'
K3=. 3 : '-(^-y)*(6+(6*y)+(3*y*y)+y^3)'
E0=. (K0 b) -K0 a
E1=. x* (K1 b) -K1 a
E2=. (*:x)* (K2 b) -K2 a
E3=. (x**:x)* (K3 b) -K3 a
EN0=. K0 xN=. ({:y)% x      NB. tau(N)-infty
ENs=. (x* K1 xN)-({:y)*EN0
E0;E1;E2;E3;EN0;ENs
)
NB.--------------------------------------------------------------------------
NB. compute emergent radiation
NB. Usage: (E;mu) IQs s;p -> I,Q,pol
IQs=: dyad define
's p'=. y
'E mu'=. x
sI=. E mm s
pI=. E mm p
m2=. *:mu
I=. sI + pI*((%3)-m2)
Q=. pI*(1 - m2)
I;Q;Q%I
)
NB.-------------------------------------------------------------------------
NB. Exponential integral function E_n(y): n EI y --> result
NB. If monadic, defaults to 1st exponential integral E_1(y)
NB. y may be any array, n should be a single value.
EIn=: (1&$:) : (dyad define)
i0=. I. 1.4>:yy=. ,y  NB. indices where  y<=1.4
i1=. I. 1.4<yy        NB. indices where  y>1.4
c0=. x&EnS i0{ yy     NB. do small arguments
c1=. x&EnCF i1{ yy    NB. do large arguments
NB. put results in original positions and restore shape:
($y)$ c1 i1} c0 i0} (#yy)#0
)
NB. For y<=1.4, we use the series expansion.
NB. n EnS y  (should be accurate to machine limit)
NB. order "n" must be single value, y may be a vector
gamma=: 0.577215664901532860606512   NB. Euler's constant
EnS=: dyad define
psi=. - gamma- +/%>: i. nm1=._1+ n=.x
b=. (psi- ^.y)* (!nm1)%~ nm1^~ -y
a=. (!k-1)* n-k=. >:i. nm1
a=. +/ a%~ (i. nm1)^~/ -y
d=. k* !kk=. n+ <: k=. >: i. 16
r=. a+ b- +/ d%~ kk^~/ -y
NB. remove E_1(0)=infinity for 1st integral:
if. n=1 do.  r=. 0 (I. y<:0)}r end.
)
NB. For y>1.4 we use a continued fraction expansion.
NB. Usage n EnCF y (error of a few times 10^-16).
NB. order "n" must be single value, y may be a vector
NB. Follows code on p 217 of Numerical Recipes
EnCF=: dyad define
nm1=. _1+ n=. x
h=. d=. %b=. n+y
c=. 1e300
i=. 0
while. (i=.>:i)<100 do.
  a=. -i*(i+nm1)
  d=. %(a*d)+ b=. b+2
  c=. b+ a%c
  h=. h*c*d
end.
h* ^-y
)
NB. ------------------------------------------------------------------------
NB. Standard memnonic definitions of J operations
mm =: +/ . *     NB. matrix multiply
DEL =: }. - }:   NB. difference between adjacent elements of vector
pwd=: 1!:43      NB. print working directory 
CD=: 1!:44       NB. change working directory
int01=: i.@>:%]  NB. int01 n --> divides [0 ... 1] into n intervals.