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 -> spline matrix operator MTsi=: 4 : 0 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 ) load '/your_path/SpMD.ijs' NB. make table of I-integrals for MTsi IIsi=: 4 : 0 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. 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.41.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 )