!     program orbitl.for

! Code converted using TO_F90 by Alan Miller
! Date: 2014-10-01  Time: 13:34:07

!     Sep. 20, 2005
!     COMPUTES MIE SCATTERING BY a modified Luneburg lens
IMPLICIT DOUBLE PRECISION (a-h,o-z)
INTEGER :: m
DIMENSION an(129),r(128),x(128),y(128)
DIMENSION ex0(128),ey0(128),fx0(128),fy0(128)
DIMENSION ex(140,128),ey(140,128),fx(140,128),fy(140,128)
DIMENSION gx(140),gg(140,128)
DIMENSION tam(128),taq(128),tad(128),tap(128)
DIMENSION tbm(128),tbq(128),tbd(128),tbp(128)
DIMENSION tem1(64),tem2(64),tem3(64),tem4(64)
DIMENSION tem5(64),tem6(64),tem7(64),tem8(64)
DIMENSION cth(361),pi(2,361)
DIMENSION srp(361),sip(361),srm(361),sim(361),pt(3,361)

pii=3.1415926536
WRITE(*,100)
100   FORMAT(' particle radius in microns=')
READ(*,101) a
a2=a*a
101   FORMAT(f14.6)
WRITE(*,102)
102   FORMAT(' magnitude of focal position/a=')
READ(*,101) f
f2=f*f
f21=f2+1.
WRITE(*,103)
103   FORMAT(' INTEGER MP <=7 FOR 2**MP interfaces=')
READ(*,104) mp
104   FORMAT(i8)
mmax=2**mp
mmax1=mmax+1
ammax=dfloat(mmax)
WRITE(*,105)
105   FORMAT(' WAVELENGTH EXAMINED IN MICRONS=')
READ(*,101) wl
WRITE(*,106)
106   FORMAT(' initial angle in degrees=')
READ(*,101) agi
WRITE(*,107)
107   FORMAT(' FINAL ANGLE IN DEGREES=')
READ(*,101) agf
nmax=361
agdel=(agf-agi)/dfloat(nmax-1)

DO  m=1,mmax
  r(m)=dfloat(m)*a/ammax
END DO
rave=.5*r(1)
an(1)=f21/f2-rave*rave/(a2*f2)
an(1)=DSQRT(an(1) )
DO  m=2,mmax
  rave=.5*(r(m)+r(m-1))
  an(m)=f21/f2-rave*rave/(a2*f2)
  an(m)=DSQRT(an(m))
END DO
an(mmax1)=1.
DO  m=1,mmax
  x(m)=2.*pii*an(m+1)*r(m) /wl
  y(m)=2.*pii*an(m)*r(m)/wl
END DO

DO  n=1,nmax
  th=pii*(agi+agdel*dfloat(n-1))/180.
  cth(n)=DCOS(th)
  srp(n)=0.
  srm(n)=0.
  sip(n)=0.
  sim(n)=0.
  pi(1,n)=0.
  pi(2,n)=1.
END DO

qlmax=2.+x(mmax)+4.315*(x(mmax)**.3333)
lmax=IFIX(SNGL(qlmax))
lstart=lmax+15

DO  m=1,mmax
  ex(lstart,m)=0.0D+00
  ey(lstart,m)=0.0D+00
  DO  ll=0,lstart-2
    l=lstart-ll
    alx=dfloat(l)/x(m)
    aly=dfloat(l)/y(m)
    ex(l-1,m)=alx-1./(ex(l,m)+ alx)
    ey(l-1,m)=aly-1./(ey(l,m)+aly)
  END DO
  alx=1./x(m)
  aly=1./y(m)
  ex0(m)=alx-1./(ex(1,m)+alx)
  ey0(m)=aly-1./(ey(1,m)+aly)
  
  fx0(m)=-DTAN(x(m))
  fy0(m)=-DTAN(y(m))
  alx=1./x(m)
  aly=1./y(m)
  fx(1,m)=-alx+1./(alx-fx0(m))
  fy(1,m)=-aly+1./(aly-fy0(m))
  DO  l=2,lmax
    alx=dfloat(l)/x(m)
    aly=dfloat(l)/y(m)
    fx(l,m)=-alx+1./(alx-fx(l-1,m))
    fy(l,m)=-aly+1./(aly-fy(l-1,m))
  END DO
END DO

DO  m=2,mmax
  gg0=DTAN(x(m-1))/DTAN(y(m))
  alx=1. /x(m-1)
  aly=1./y(m)
  gg(1,m)=gg0*((alx+fx(1,m-1))/(aly+fy(1,m)))
  gg(1,m)=gg(1,m)*((aly+ey(1,m))/(alx + ex(1,m-1)))
  DO  l=2,lmax
    alx=dfloat(l)/x(m-1)
    aly=dfloat(l)/y(m)
    gg(l,m)=gg(l-1,m)*((alx+fx(l,m-1))/(aly+fy(l,m)))
    gg(l,m)=gg(l,m)*((aly+ey(l,m) ) / (alx + ex(l,m-1)))
  END DO
END DO

gx0=-DTAN(x(mmax))
alx=1./x (mmax)
gx(1) =gx0*((alx+fx(1,mmax))/(alx+ex(1,mmax)))
DO  l=2,lmax
  alx=dfloat(l) /x(mmax)
  gx(l)=gx(l-1)*((alx+fx(l,mmax))/(alx + ex(l,mmax)))
END DO

DO  l=1,lmax
  q=dfloat(l)
  q1=q+1.
  q21=2.*q+1.
  qq1=q21/q
  qq2=q1/q
  qq3=q21/(q*q1)
  
  DO  m=1,mmax
    tam(m)=an(m+1)*ey(l,m)-an(m) *ex(l,m)
    taq(m)=an(m+1)*fy(l,m)-an(m) *fx(l,m)
    tad(m)=an(m+1)*ey(l,m)-an(m) *fx(l,m)
    tap(m)=an(m+1)*fy(l,m)-an(m) *ex(l,m)
    tbm(m)=an(m)*ey(l,m) -an(m+1)*ex(l,m)
    tbq(m)=an(m)*fy(l,m) -an(m+1)*fx(l,m)
    tbd(m)=an(m)*ey(l,m) -an(m+1)*fx(l,m)
    tbp(m)=an(m)*fy(l,m) -an(m+1)*ex(l,m)
    
  END DO
  
  DO  j=1,mp
    jp2=2**(mp-j)
    ij1 = 2** (j-1)
    
    DO  k=1,jp2
      k2=2*k
      k21=k2-1
      kk=ij1*k21+1
      grat=gg(l,kk)
      tem1(k)=tad(k21)*tam(k2)-tam(k21)*tap(k2)*grat
      tem2(k)=taq(k21)*tad(k2)-tap(k21)*taq(k2)*grat
      tem3(k)=tad(k21)*tad(k2)-tam(k21)*taq(k2)*grat
      tem4(k)=taq(k21)*tam(k2)-tap(k21)*tap(k2)*grat
      tem5(k)=tbd(k21)*tbm(k2)-tbm(k21)*tbp(k2)*grat
      tem6(k)=tbq(k21)*tbd(k2)-tbp(k21)*tbq(k2)*grat
      tem7(k)=tbd(k21)*tbd(k2)-tbm(k21)*tbq(k2)*grat
      tem8(k)=tbq(k21)*tbm(k2)-tbp(k21)*tbp(k2)*grat
    END DO
    
    DO  k=1,jp2
      tam(k)=tem1(k)
      taq(k)=tem2(k)
      tad(k)=tem3(k)
      tap(k)=tem4(k)
      tbm(k)=tem5(k)
      tbq(k)=tem6(k)
      tbd(k)=tem7(k)
      tbp(k)=tem8(k)
    END DO
  END DO
  
  tam(1)=tam(1)*gx(l)
  tbm(1)=tbm(1)*gx(l)
  aden=tam(1) *tam (1) +tad (1) *tad (1)
  ar=tam(1)*tam(1)/aden
  ai=-tam(1)*tad(1)/aden
  bden=tbm(1)*tbm(1)+tbd(1)*tbd(1)
  br=tbm(1)*tbm(1)/bden
  bi=-tbm(1)*tbd(1)/bden
  abrp=qq3*(ar+br)
  abip=qq3*(ai+bi)
  abrm=qq3*(ar-br)
  abim=qq3*(ai-bi)
  
  DO  n=1,nmax
    ta=q*cth(n)*pi(2,n)-q1*pi(1,n)
    pitap=pi(2,n)+ta
    pitam=pi(2,n)-ta
    srp(n)=srp(n)+abrp*pitap
    srm(n)=srm(n)+abrm*pitam
    sip(n)=sip(n)+abip*pitap
    sim(n)=sim(n)+abim*pitam
    temp=pi(2,n)
    pi(2,n)=qq1*cth(n)*pi(2,n)-qq2*pi(1,n)
    pi(1,n)=temp
  END DO
END DO

DO  n=1,nmax
  s1r=.5*(srp(n)+srm(n))
  s2r=.5*(srp(n)-srm(n))
  s1i=.5*(sip(n)+sim(n))
  s2i=.5*(sip(n)-sim(n))
  pt(1,n)=s1r*s1r+s1i*s1i
  pt(2,n)=s2r*s2r + s2i*s2i
  pt(3,n)= .5*(pt(1,n)+pt(2,n))
END DO

OPEN(UNIT=1,FILE='OUT.DAT')

DO  n=1,nmax
  th=agi+agdel*dfloat(n-1)
  
  WRITE(1,998) th,pt(1,n),pt(2,n),pt(3,n)
  998   FORMAT(f14.1,3E15.5)
END DO
CLOSE(UNIT=1,STATUS='KEEP')
400   CONTINUE
END
