; June 18, 2002
;
; This list of IDL routines contains two main routines which
; implement the equations in Agol (2002):
;
;  a) limbdark - computes the lightcurve of a microlensing/
;     occultation event including limb-darkening of the background
;     star
;
;  b) extsource- computes the lightcurve of a microlensing/occultation
;     event for a uniform source
;
;  The input and output parameters for each routine are described
;  at the start of each routine.
;
;  There are several other routines in this file which are called
;  by the above three routines.  The elliptical integral routines
;  are versions of the routines in Numerical Recipes in FORTRAN
;  by Press, Teukolsky, Vetterling, and Flannery (1992) modified for IDL.
;
;  To run in unix or linux:
;     1)  Save this file as, e.g., 'microccult.pro'
;     2)  Start IDL
;     3)  Start an x-window:  window,0,retain=2
;     3)  From the IDL prompt, type '.run microccult', twice.
;         You should see a list of compiled routines.
;     4)  Run one of the routines.
;
;  Example calls:
;
;  a) First, an example of using limbdark:
;
;    limbdark,1.d0,1.d0,0.4d0,0.4d0,b0,mulimb,1
;
;  You should see the following lines:
;  IDL>  limbdark,1.d0,1.d0,0.4d0,0.4d0,b0,mulimb,1
;  rs       1.00000 rl       1.00000
;  rl > 1, occulted
;  % Program caused arithmetic error: Floating divide by 0
;  % Program caused arithmetic error: Floating illegal operand
;
;  and a plot should appear in the X-window showing a dashed line
;  for a uniform, extended source, a solid line for a limb-darkened
;  source, and a diamond for the analytic limb-darkened amplification
;  for a source on axis.
; 
;  If you would like to use different impact parameters, then
;  simply create a vector b0 and then call limbdark.  The limb-darkened
;  amplification will be output as mulimb.  To change the plot axes,
;  then simply add IDL plotting parameters at the end of the call.
;  
;  b)  Next, an example of using extsource:
;
;  IDL> extsource,1.d0,0.95d0,1,1,1,b0,muom,muop,mu,/noerase,ys=1,yr=[0.,3.]
;  % Program caused arithmetic error: Floating divide by 0
;
;  You should see a plot with three lines:  solid line is for extended
;  microlensed/occulted source, dotted line is for a point source
;  while dashed is for a microlensed extended source.
; 
;  If you make use of these routines, please cite Agol (2002) and please
;  drop me a note at agol AT astro.washington.edu.  If you have any problems
;  with the routines, please let me know.
;  
;  Here are the routines:

pro limbdark,rs,rl,u1,u2,b0,mulimb,plotquery,_extra=e
; This routine uses the results for a uniform source to
; compute the lightcurve for a limb-darkened source
; (5-1-02 notes)
;Input:
;  rs        radius of the source in units of the Einstein radius
;  rl        radius of the lens   in units of the Einstein radius
;  u1        linear limb-darkening coefficient
;  u2        quadratic limb-darkening coefficient
;  b0        impact parameter normalized to source radius
;  plotquery =1 to plot magnification,  =0 for no plots
;  _extra=e  plot parameters
;Output:
; mulimb limb-darkened magnification
; 
; First, make grid in radius:
; Call magnification of uniform source:
extsource,rs,rl,0,0,0,bvec=b0,bt0,muom,muop,mu
mulimb0=muom+muop
mulimb=mulimb0
fac=max(abs(mulimb0-1.))
;print,'fac ',fac
print,'rs ',rs,' rl ',rl
omega=1.-u1/3.-u2/6.
alp1=2.*(u1+2.*u2)
alp2=3.*(1.-u1-u2)*rs-3.*u2/2./rs*(2.+rs^2)
nb=n_elements(b0)
for j=0,nb-1 do begin
nr=2
dmumax=1.
while (dmumax gt fac*4.d-3) do begin
  mulimb1=mulimb(j)
  nr=nr*2
  dt=0.5*!pi/double(nr)
  t=dt*dindgen(nr+1)
  th=t+0.5*dt
  r=rs*sin(t)
  mulimb(j)=cos(th(nr-1))^2*rs^2*mulimb0(j)/(rs-r(nr-1))
  mulimb2=cos(th(nr-1))^3*rs^2*mulimb0(j)/(rs-r(nr-1))
  for i=1,nr-1 do begin
; Calculate uniform magnification at intermediate radii:
    extsource,r(i),rl,0,0,0,bvec=b0(j)*rs/r(i),bt,muom,muop,mu
    mu=muom+muop
; Equation (26):
    mulimb(j)= mulimb(j) +r(i)^2*mu*(cos(th(i-1))^2/(r(i)-r(i-1))-cos(th(i))^2/(r(i+1)-r(i)))
    mulimb2=mulimb2+r(i)^2*mu*(cos(th(i-1))^3/(r(i)-r(i-1))-cos(th(i))^3/(r(i+1)-r(i)))
  endfor
  mulimb(j)=((1.-u1-u2)*mulimb0(j)+mulimb(j)*(u1+2.*u2)*dt/rs-mulimb2*u2*dt/rs)/omega
  indx=where((mulimb(j) ne 0.) and (mulimb1 ne 0.))
  if(mulimb(j) ne 0. and mulimb1 ne 0.) then begin
    dmumax=abs(mulimb(j)-mulimb1)/(mulimb(j)+mulimb1)
  endif else begin
    dmumax=0.
  endelse
;  print,'Difference ',dmumax,' nr ',nr
endwhile
endfor
if(plotquery eq 1) then plot,bt0,mulimb,_extra=e
if(plotquery eq 1) then oplot,bt0,mulimb0,linestyle=2
; Plot the analytic function for zeta_0=0 (5-9-2002 notes):
k=rs/sqrt(4.+rs^2)
; 
if(rl eq 0. or rl*rs lt 1.-rl^2) then begin
  print,'unocculted'
; Equation (29)
  mu0=(alp1/3./k*((2.+rs^2)*elle(0.5*!pi,k)-2.*ellf(0.5*!pi,k))+alp2*rs/3./k+$
       4.*u2/rs^2*asinh(rs/2.))/omega/rs^2
endif
if(rl lt 1. and rl*rs ge 1.-rl^2) then begin
  print,'rl <1, occulted'
  bl=1./rl-rl
  alp3=rs^2-bl^2
  phi3=acos(-bl/rs)
  ee=2.*elle(0.5*!pi,k)-elle(phi3,k)
  ff=2.*ellf(0.5*!pi,k)-ellf(phi3,k)
  s=(1.-rl)/abs(1.-rl)
; Equation (27)
  mu0=sqrt(4.+rs^2)/6./omega/rs^3*(alp1*((2.+rs^2)*ee-2.*ff)+alp2*rs)+$
       (bl*sqrt(4.+bl^2)+alp3)*(alp1*sqrt(alp3)+alp2-3.*u2/2./rs*alp3)/6./omega/rs^3+$
       2.*u2/rs^4/omega*(alp3/8.*(2.+rs^2)+asinh(rs/2.)+asinh(bl/2.))
endif
if(rl ge 1. and rs*rl ge rl^2-1.) then begin
  print,'rl > 1, occulted'
  bl=rl-1./rl
  alp3=rs^2-bl^2
  phi3=acos(bl/rs)
  ee=elle(phi3,k)
  ff=ellf(phi3,k)
; Equation (27)
  mu0=sqrt(4.+rs^2)/6./omega/rs^3*(alp1*((2.+rs^2)*ee-2.*ff)+alp2*rs)+$
       (-bl*sqrt(4.+bl^2)+alp3)*(alp1*sqrt(alp3)+alp2-3.*u2/2./rs*alp3)/6./omega/rs^3+$
       2.*u2/rs^4/omega*(alp3/8.*(2.+rs^2)+asinh(rs/2.)+asinh(-bl/2.))
endif
if(rl ge 1. and rs*rl le rl^2-1.) then mu0=0.
;print,mu0,mulimb(0)
if(plotquery eq 1) then oplot,[0.,0.],[mu0,mu0],psym=4,symsize=1,/noclip
b0=bt0
return
end

function asinh,x
return,alog(x+sqrt(x^2+1.))
end

pro extsource,rs,rl,pc,ec,oc,bmin=bmin,bmax=bmax,bvec=bvec,b0,muom,muop,mu,_extra=e
; This routine computes the lightcurve for microlensing & occultation of a uniform
; extended source by a point lens (Agol 2002).
;
;Input:
;  rs      source size in units of Einstein radius projected onto source plane
;  rl      radius of the lens in units of the Einstein radius
;  bvec    position of source in source-radius units
;  bmin    minimum impact parameter (if bvec is unspecified)
;  bmax    maximum impact parameter (if bvec is unspecified)
;  ec      =1 for plotting extended (unocculted) magnification =0 no plot
;  pc      =1 for plotting point-source,occulted magnification =0 no plot
;  oc      =1 for plotting extended,occulted magnification =0 no plot
;
;Output:
;  b0      impact parameter grid
;  muom    magnification of negative (inner) image
;  muop    magnification of positive (outer) image
;  mu      1/2 of the magnification of an unocculted extended uniform source
; 
if(n_elements(bmin) eq 0) then bmin=0.d0
if(n_elements(bmax) eq 0) then bmax=4.d0
; If the call doesn't specify a list of impact parameters, then 
; generate a list:
if(n_elements(bvec) eq 0) then begin
  nb=200
  bvec=(bmax-bmin)*dindgen(nb+1)/double(nb)+bmin
endif else begin
  nb=n_elements(bvec)-1
endelse
; position of source in Einstein-radius units (in the paper
; this is the parameter \zeta_0):
be=bvec*rs
mu=dblarr(n_elements(be))
; Now, compute expression from Witt & Mao (1994):
n=4.d0*be*rs/(be+rs)^2
k=sqrt(4.d0*n/(4.d0+(be-rs)^2))
indx=where(abs(be-rs) gt 1.d-7)
mu(indx)=gofphi(!pi/2.,be(indx),rs)/4./!pi/rs^2
indx=where(abs(be-rs) lt 1.d-7)
if(indx gt 0.d0) then begin
; Use special formula when the impact parameter equals the source star radius:
  mu(indx)=(rs+(1.d0+rs^2)*atan(rs))/!pi/rs^2
endif
; plot total extended-source (unocculted) magnification:
if(ec eq 1) then plot,bvec,2.d0*mu,linestyle=2,_extra=e

; Next, plot point-source expression:
if(rl lt 1.d0) then begin
  mup=0.5d0*( 1.d0+(be^2+2.d0)/be/sqrt(be^2+4.d0))
  bl=1.d0/rl-rl
  mum=0.5d0*(-1.d0+(be^2+2.d0)/be/sqrt(be^2+4.d0))*(be le bl)
endif else begin
  bl=rl-1.d0/rl
  mup=0.5d0*( 1.d0+(be^2+2.d0)/be/sqrt(be^2+4.d0))*(be gt bl)
  mum=0.d0*bvec
endelse
if(pc eq 1) then plot,bvec,mup+mum,linestyle=1,_extra=e

; Now, compute the occulted magnification:
; First, set the magnifications equal to the unocculted magnification:
muom=mu-0.5d0
muop=mu+0.5d0
; There are 3 cases - rl=0, 0<rl<1, rl>1:
; (1) rl = 0
; There is no occultation - same as Witt & Mao (1995) expression:
if(rl eq 0.d0) then  goto,finish

if(rl eq 1.d0) then begin
; Inner image is occulted, outer is not:
  muom=0.*muom
  goto,finish
endif
; (2) 0 < rl < 1:
;     In this case, the inner (negative parity) image can be occulted, 
;     while the outer (positive parity) image is never occulted:
if(rl lt 1.d0) then begin
; Project lens radius onto source plane:
  bl=1.d0/rl-rl
for i=0,nb do begin
; (a) Inner image is completely occulted by lens star:
  if(be(i) gt bl+rs) then begin
    muom(i)=0.d0 
    goto,positive
  endif
; (b) When the edge of the star crosses the lens singularity, we need a 
; special expression:
  if(abs(be(i)-rs) lt 1.d-7) then begin
; Partial occultation (equation 19):
    if(rs gt 0.5d0*bl) then begin
      v1=sqrt(4.d0+bl^2)
      v2=sqrt(4.d0*rs^2-bl^2)
      a1=(0.25d0*v2*(bl-v1)+(1.d0+rs^2)*(acos(0.5d0*bl/rs)-atan(v2/v1)))/!pi/rs^2
      a2=rl^2/!pi/rs^2*acos(0.5d0*bl/rs)
      muom(i)=muom(i)+a1-a2
    endif
    goto,positive
  endif
; (c) Inner image is unocculted - use Witt & Mao expression:
  if(be(i) le bl-rs) then goto,positive
; (d) When the source star covers the unocculted window, then the magnification
; is a constant:
  if(be(i) lt rs-bl) then begin
    muom(i)=(1.d0-rl^2)/rs^2 
    goto,positive
  endif
;  if(abs(bl-rs) lt be(i) and be(i) lt bl+rs) then begin
; (e) The inner image is partially occulted by the lens star:
; (see 4/18/02 notes) equation (16)
    u0=bl^2
    u1=(be(i)-rs)^2
    u2=(be(i)+rs)^2
    u3=be(i)^2-rs^2
; Now, use notation from equation (18) in paper:
    phi0=acos(sqrt(u1*(u2-u0)/u0/(u2-u1)))
    phi1=acos((u1+u2-2.*u0)/(u2-u1))
    phi2=acos((u3+u0)/2./be(i)/bl)
    muom(i)=1.d0/rs^2*(be(i) lt rs)-(-4.*u3/abs(u3)*phi0+2.*(1.+rs^2)*phi1+$
       4.*rl^2*phi2+sqrt((u2-u0)*(u0-u1))*(sqrt(1.+4./u0)-1.)-$
       gofphi(phi0,be(i),rs))/4./!pi/rs^2
positive:
; The outer (positive image) is unocculted.
  muop(i)=muop(i)
endfor
goto,finish
endif 

; Now for rl gt 1 case:
; Compute the size of the lens star projected onto the source plane:
bl=rl-1.d0/rl
for i=0,nb do begin
; Image is unocculted - use usual formula:
  if(be(i) gt bl+rs) then goto, negative
; Image is completely occulted:
  if(rs le bl and be(i) le bl-rs) then begin
    muop(i)=0.d0
    goto,negative
  endif
; Image has a hole in the center:
  if(rs gt bl and abs(be(i)) lt rs-bl) then begin
    muop(i)=mu(i)+0.5d0-(rl^2-1.d0)/rs^2 
    goto,negative
  endif
; Image is partially occulted:  use equations (22) & (23):
  u1=(be(i)-rs)^2
  u2=(be(i)+rs)^2
  u3=be(i)^2-rs^2
  u0=bl^2
  psi0=acos(sqrt((u0-u1)*(4.+u2)/(4.+u0)/(u2-u1)))
  phi0=acos(sqrt(u1*(u2-u0)/u0/(u2-u1)))
  phi1=acos((u1+u2-2.*u0)/(u2-u1))
  phi2=acos((u3+u0)/2./be(i)/bl)
  psi1=!pi/2.-phi0
  psi2=!pi-phi1+2.*acos(sqrt(u0*(4.+u0)/(u0*(4+u1+u2)-u1*u2)))
  if(u3 ne 0.d0) then begin
    muop(i)=(-4.*u3/abs(u3)*psi1+2.*(1.+rs^2)*psi2-4.*rl^2*phi2+$
        sqrt((u2-u0)*(u0-u1))*(sqrt(u0/(4.+u0))+1.)+gofphi(psi0,be(i),rs))/4./!pi/rs^2
  endif else begin
    muop(i)=(-4.*psi1+2.*(1.+rs^2)*psi2-4.*rl^2*phi2+$
        sqrt((u2-u0)*(u0-u1))*(sqrt(u0/(4.+u0))+1.)+gofphi(psi0,be(i),rs))/4./!pi/rs^2
  endelse
negative:
; Negative (inner) image is always occulted:
  muom(i)=0.d0
endfor
finish:
if(oc eq 1) then plot,bvec,muop+muom,linestyle=0,_extra=e
b0=bvec
return
end

function gofphi,phi,z0,rs
u1=(z0-rs)^2
u2=(z0+rs)^2
u3=z0^2-rs^2
en=1.-u1/u2
k=sqrt(4.*(u2-u1)/u2/(4.+u1))
;ee=elle(phi,k)
;ff=ellf(phi,k)
;pp=ellpi(phi,-en,k)
one=1.d0+dblarr(n_elements(z0))
sp=sin(phi)*one
cp2=1.d0-sp*sp
yy=1.d0-k*k*sp*sp
ff=sp*rf(cp2,yy,one)
ee=ff-k*k*sp*sp*sp/3.d0*rd(cp2,yy,one)
pp=ff+en*sp*sp*sp/3.d0*rj(cp2,yy,one,1.d0-en*sp*sp)
gf=(u2*(4.+u1)*ee-(u1*u2+8.*u3)*ff+4.*u1*(1.+rs^2)*pp)/sqrt(u2*(4.+u1))
return,gf
end

FUNCTION elle,phi,ak
;REAL elle,ak,phi
;REAL cc,q,s,rd,rf
if(n_elements(phi) eq 1) then phi=phi+fltarr(n_elements(ak))
s=sin(phi)
cc=cos(phi)^2
q=(1.-s*ak)*(1.+s*ak)
one=1.+fltarr(n_elements(ak))
ellee=s*(rf(cc,q,one)-((s*ak)^2)*rd(cc,q,one)/3.)
return,ellee
END

function ellf,phi,ak
;REAL ellf,ak,phi
;USES rf
;REAL s,rf
if(n_elements(phi) eq 1) then phi=phi+fltarr(n_elements(ak))
one=1.+fltarr(n_elements(ak))
s=sin(phi)
ellf=s*rf(cos(phi)^2,(1.-s*ak)*(1.+s*ak),one)
return,ellf
END
FUNCTION ellpi,phi,en,ak
;REAL ellpi,ak,en,phi
;USES rf,rj
;REAL cc,enss,q,s,rf,rj
if(n_elements(phi) eq 1) then phi=phi+dblarr(n_elements(ak))
s=sin(phi)
enss=en*s*s
cc=cos(phi)^2
q=(1.d0-s*ak)*(1.d0+s*ak)
one=1.d0+dblarr(n_elements(ak))
rjj=rj(cc,q,one,one+enss)
ellpi=s*(rf(cc,q,one)-enss*rjj/3.d0)
return,ellpi
END

function rc,x,y
;REAL rc,x,y,ERRTOL,TINY,SQRTNY,BIG,TNBG,COMP1,COMP2,THIRD,C1,C2,
;*C3,C4
ERRTOL=.04
;ERRTOL=1.d-3
TINY=1.69e-38
SQRTNY=1.3e-19
BIG=3.E37
TNBG=TINY*BIG
COMP1=2.236/SQRTNY
COMP2=TNBG*TNBG/25.
THIRD=1./3.
C1=.3
C2=1./7.
C3=.375
C4=9./22.
;REAL alamb,ave,s,w,xt,yt
rcc=dblarr(n_elements(x))
if(n_elements(x) ne n_elements(y)) then begin
   print,'wrong size in rc'
   return,dblarr(n_elements(x))
endif
for i=0,n_elements(x)-1 do begin
  if(x(i) lt 0. or y(i) eq 0. or (x(i)+abs(y(i))) lt TINY $
     or (x(i)*abs(y(i))) gt BIG or (y(i) lt -COMP1 and x(i) $
     gt 0. and x(i) lt COMP2)) then begin
      print,'invalid arguments in rc'
;      return,dblarr(n_elements(x))
     rcc(i)=0.d0
  endif else begin
;if(x.lt.0..or.y.eq.0..or.(x+abs(y)).lt.TINY.or.(x+
;     *abs(y)).gt.BIG.or.(y.lt.-COMP1.and.x.gt.0..and.x.lt.COMP2))pause 
;     *'invalid arguments in rc'
  if(y(i) gt 0.)then begin
    xt=x(i)
    yt=y(i)
    w=1.
  endif else begin
    xt=x(i)-y(i)
    yt=-y(i)
    w=sqrt(x(i))/sqrt(xt)
  endelse
LAB1:
  alamb=2.*sqrt(xt)*sqrt(yt)+yt
  xt=.25*(xt+alamb)
  yt=.25*(yt+alamb)
  ave=THIRD*(xt+yt+yt)
  s=(yt-ave)/ave
  if(abs(s) gt ERRTOL) then goto,LAB1
  rcc(i)=w*(1.+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave)
  endelse
endfor
return,rcc
END

function rd,x,y,z
;      REAL rd,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6
;      PARAMETER (ERRTOL=.05,TINY=1.e-25,BIG=4.5E21,C1=3./14.,C2=1./6.,
;     *C3=9./22.,C4=3./26.,C5=.25*C3,C6=1.5*C4)
;      REAL alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty,
;     *sqrtz,sum,xt,yt,zt
ERRTOL=.05d0
;ERRTOL=1.d-3
TINY=1.d-25
BIG=4.5d21
C1=3.d0/14.d0
C2=1.d0/6.d0
C3=9.d0/22.d0
C4=3.d0/26.d0
C5=.25d0*C3
C6=1.5d0*C4
rdd=dblarr(n_elements(x))
if(n_elements(x) ne n_elements(y) or n_elements(x) ne n_elements(z)) then begin
   print,'wrong size in rd'
   return,dblarr(n_elements(x))
endif
for i=0,n_elements(x)-1 do begin
  if(min([x(i),y(i),z(i)]) lt 0.d0 or min([x(i)+y(i),x(i)+z(i),y(i)+z(i)]) lt TINY $
     or max([x(i),y(i),z(i)]) gt BIG) then begin
      print,'invalid arguments in rd'
;      return,dblarr(n_elements(x))
     rdd(i)=0.d0
  endif else begin
  xt=x(i)
  yt=y(i)
  zt=z(i)
  sum=0.d0
  fac=1.d0
LAB1:
  sqrtx=sqrt(xt)
  sqrty=sqrt(yt)
  sqrtz=sqrt(zt)
  alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
  sum=sum+fac/(sqrtz*(zt+alamb))
  fac=.25d0*fac
  xt=.25d0*(xt+alamb)
  yt=.25d0*(yt+alamb)
  zt=.25d0*(zt+alamb)
  ave=.2d0*(xt+yt+3.d0*zt)
  if(ave eq 0.d0) then begin
    delx=0.d0
    dely=0.d0
    delz=0.d0
  endif else begin
    delx=(ave-xt)/ave
    dely=(ave-yt)/ave
    delz=(ave-zt)/ave
  endelse
  if(max([abs(delx),abs(dely),abs(delz)]) gt ERRTOL) then goto,LAB1
  ea=delx*dely
  eb=delz*delz
  ec=ea-eb
  ed=ea-6.d0*eb
  ee=ed+ec+ec
  rdd(i)=3.d0*sum+fac*(1.d0+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3*$
      ec+delz*C4*ea)))/(ave*sqrt(ave))
  endelse
endfor
return,rdd
END

; Program for evaluating elliptic integrals: in notation of
; Gradsteyn and Rhyzik.
function rf,x,y,z
;double precision FUNCTION rf(x,y,z)
;implicit none
;double precision rf,x,y,z,ERRTOL,TINY,BIG,THIRD,C1,C2,C3,C4
ERRTOL=0.0025d0
;ERRTOL=1.d-5
TINY=1.5d-38
BIG=3.d37
THIRD=1.d0/3.d0
C1=1.d0/24.d0
C2=0.1d0
C3=3.d0/44.d0
C4=1.d0/14.d0
rff=dblarr(n_elements(x))
if(n_elements(x) ne n_elements(y) or n_elements(x) ne n_elements(z)) then begin
   print,'wrong size in rf'
   return,dblarr(n_elements(x))
endif
for i=0,n_elements(x)-1 do begin
  if(min([x(i),y(i),z(i)]) lt 0.d0 or min([x(i)+y(i),x(i)+z(i),y(i)+z(i)]) lt TINY $
     or max([x(i),y(i),z(i)]) gt BIG) then begin
      print,'invalid arguments in rf'
      print,min([x(i),y(i),z(i)]),min([x(i)+y(i),x(i)+z(i),y(i)+z(i)]),max([x(i),y(i),z(i)])
;      return,dblarr(n_elements(x))
      rff(i)=0.d0
  endif else begin
  xt=x(i)
  yt=y(i)
  zt=z(i)
LAB1:
  sqrtx=sqrt(xt)
  sqrty=sqrt(yt)
  sqrtz=sqrt(zt)
  alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
  xt=0.25d0*(xt+alamb)
  yt=0.25d0*(yt+alamb)
  zt=0.25d0*(zt+alamb)
  ave=THIRD*(xt+yt+zt)
  if(ave eq 0.d0) then begin
    delx=0.d0
    dely=0.d0
    delz=0.d0
  endif else begin
    delx=(ave-xt)/ave
    dely=(ave-yt)/ave
    delz=(ave-zt)/ave
  endelse
  if(max([abs(delx),abs(dely),abs(delz)]) gt ERRTOL) then goto,LAB1
  e2=delx*dely-delz*delz
  e3=delx*dely*delz
  rff(i)=(1.d0+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave)
  endelse
endfor
return,rff
END

function rj,x,y,z,p
;USES rc,rf
;REAL a,alamb,alpha,ave,b,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee,
;*fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt,rc,rf
;REAL rj,p,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6,C7,C8
ERRTOL=.05
;ERRTOL=1.d-3
TINY=2.5e-13
BIG=9.E11
C1=3./14.
C2=1./3.
C3=3./22.
C4=3./26.
C5=.75*C3
C6=1.5*C4
C7=.5*C2
C8=C3+C3
rjj=dblarr(n_elements(x))
if(n_elements(x) ne n_elements(y) or n_elements(x) ne n_elements(z) or $
      n_elements(x) ne n_elements(p)) then begin
   print,'wrong size in rj'
   return,dblarr(n_elements(x))
endif
for i=0,n_elements(x)-1 do begin
  if(min([x(i),y(i),z(i)]) lt 0.d0 or min([x(i)+y(i),x(i)+z(i),y(i)+z(i),abs(p(i))]) lt TINY $
     or max([x(i),y(i),z(i),abs(p(i))]) gt BIG) then begin
      rjj(i) = 0.d0
      print,'invalid arguments in rj'
      print, x(i),y(i),z(i),p(i)
;      return,dblarr(n_elements(x))
  endif else begin
  sum=0.
  fac=1.
  if(p(i) gt 0 )then begin
    xt=x(i)
    yt=y(i)
    zt=z(i)
    pt=p(i)
  endif else begin
    xt=min([x(i),y(i),z(i)])
    zt=max([x(i),y(i),z(i)])
    yt=x(i)+y(i)+z(i)-xt-zt
    a=1./(yt-p(i))
    b=a*(zt-yt)*(yt-xt)
    pt=yt+b
    rho=xt*zt/yt
    tau=p(i)*pt/yt
    rcx=rc(rho,tau)
  endelse
LAB1:
  sqrtx=sqrt(xt)
  sqrty=sqrt(yt)
  sqrtz=sqrt(zt)
  alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
  alpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)^2
  beta=pt*(pt+alamb)^2
  sum=sum+fac*rc(alpha,beta)
  fac=.25*fac
  xt =.25*(xt+alamb)
  yt =.25*(yt+alamb)
  zt =.25*(zt+alamb)
  pt =.25*(pt+alamb)
  ave=.2*(xt+yt+zt+pt+pt)
  delx=(ave-xt)/ave
  dely=(ave-yt)/ave
  delz=(ave-zt)/ave
  delp=(ave-pt)/ave
  if(max([abs(delx),abs(dely),abs(delz),abs(delp)]) gt ERRTOL)then goto,LAB1
  ea=delx*(dely+delz)+dely*delz
  eb=delx*dely*delz
  ec=delp^2
  ed=ea-3.*ec
  ee=eb+2.*delp*(ea-ec)
  rjj(i)=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))$
      +delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave))
  if (p(i) le 0.) then rjj(i)=a*(b*rjj(i)+3.*(rcx-rf(xt,yt,zt)))
  endelse
endfor
return,rjj
END