C     program orbitl.for
C     Sep. 20, 2005
C     COMPUTES MIE SCATTERING BY a modified Luneburg lens
C     IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      IMPLICIT REAl(KIND=8)  (A-H,O-Z)


      integer m
      DIMENSION AN(513),r(512),X(512),Y(512)
      DIMENSION ex0(512),ey0(512),fx0(512),fy0(512)
CC    lmax = 800 with r = 60 microns
      DIMENSION ex(800,512),ey(800,512),fx(800,512),fy(800,512)
      DIMENSION gx(800),gg(800,512)
      DIMENSION TAM(512),TAQ(512),TAD(512),TAP(512)
      DIMENSION TBM(512),TBQ(512),TBD(512),TBP(512)
      DIMENSION TEM1(256),TEM2(256),TEM3(256),TEM4(256)
      DIMENSION TEM5(256),TEM6(256),TEM7(256),TEM8(256)
CC      DIMENSION CTH(361),PI(2,361)
CC      DIMENSION SRP(361),SIP(361),SRM(361),SIM(361),PT(3,361)
      DIMENSION CTH(1801),PI(2,1801)
      DIMENSION SRP(1801),SIP(1801),SRM(1801),SIM(1801),PT(3,1801)
C
      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 <=9 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
      
CC    thw
CC      real :: start, finish
      call cpu_time(start)
CC    ! put code to test here
CC      NMAX=361
      NMAX=1801
      AGDEL=(AGF-AGI)/DFLOAT(NMAX-1)
C
        DO 200 M=1,MMAX
      r(m)=dfloat(m)*a/ammax
200   CONTINUE
      rave=.5*r(1)
      an(1)=f21/f2-rave*rave/(a2*f2)
      an(1)=dsqrt(an(1) )
        do 201 m=2,mmax
      rave=.5*(r(m)+r(m-1))
      an(m)=f21/f2-rave*rave/(a2*f2)
      an(m)=dsqrt(an(m))
201   continue
      AN(MMAX1)=1.
        DO 202 M=1,MMAX
      X(M)=2.*PII*AN(M+1)*R(m) /WL
      Y(M)=2.*PII*AN(M)*R(m)/WL
202   CONTINUE
C
        do 211 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.
211   continue
C
      QLMAX=2.+x(MMAX)+4.315*(x(MMAX)**.3333)
      LMAX=ifix(SNGL(QLMAX))
      lstart=lmax+15
C
        DO 203 M=1,MMAX
      ex(lstart,m)=0.0d+00
      ey(lstart,m)=0.0d+00
        do 215 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)
215   continue
      alx=1./x(m)
      aly=1./y(m)
      ex0(m)=alx-1./(ex(1,m)+alx)
      ey0(m)=aly-1./(ey(1,m)+aly)
C
      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 216 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))
216   continue
203   CONTINUE
C
        do 219 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 217 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)))
217   continue
219   continue
c
      gx0=-dtan(x(mmax))
      alx=1./x (mmax)
      gx(1) =gx0*((alx+fx(1,mmax))/(alx+ex(1,mmax)))
        do 218 L=2,lmax
      alx=dfloat(L) /x(mmax)
      gx(L)=gx(L-1)*((alx+fx(L,mmax))/(alx + ex(L,mmax)))
218   continue
C
        do 204 L=1,LMAX
      q=Dfloat(L)
      Q1=Q+1.
      q21=2.*q+1.
      QQ1=Q21/Q
      QQ2=Q1/Q
      QQ3=Q21/(Q*Q1)
C
        DO 205 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)

205   CONTINUE
C
        DO 206 J=1,MP
      JP2=2**(MP-J)
      ij1 = 2** (j-1)
C
        DO 207 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
207   CONTINUE
C
        DO 208 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)
208   CONTINUE
206   CONTINUE
C
      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)
C
        DO 212 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
212   CONTINUE
204   CONTINUE
C
        do 213 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))
213   continue
C
      OPEN(UNIT=1,FILE='OUT.DAT')
C     thw
      anorm=PT(1,1)
      
        DO 210 N=1,NMAX
      TH=AGI+AGDEL*DFLOAT(N-1)
C     thw
      PT(1,N)=PT(1,N)/anorm*255.
      PT(2,N)=PT(2,N)/anorm*255.
      PT(3,N)=PT(3,N)/anorm*255.
      
      WRITE(1,998) TH,PT(1,N),PT(2,N),PT(3,N)
998   FORMAT(F14.1,3E18.8)
210   CONTINUE
      CLOSE(UNIT=1,STATUS='KEEP')
400   continue
CC    thw
      call cpu_time(finish)
      print '("Time = ",f6.3," seconds.")',finish-start
      print '("lstart = ",I6)',lstart
      print '("mmax = ",I6)',mmax
      print '("lmax = ",I6)',lmax
      read(*,101) aaaaa
      end
