! ZoneSize3.zpl
!
! Written by Dan Brown, MEMS Optical, Huntsville, AL 35806
! 3/18/99
! Modified 4/7/99
!
! Computes and prints diffractive zone width across binary 1 or 2
! Scan distance is set by semidiameter in lens data editor.
!-----------------------------------------------------------------
file$ = $FILEPATH()
PRINT file$
IF UNIT()==0
   const = 6283.185
   un$ = "(mm)"
ELSE
   IF UNIT()==1
      const = 62831.85
      un$ = "(cm)"
   ELSE
      IF UNIT()==2
         const = 159592.9
         un$ = "(in)"
      ELSE
         IF UNIT()==3
            const = 6283185
            un$ = "(m)"
         ENDIF
      ENDIF
   ENDIF
ENDIF
INPUT "Enter surface number of Binary 1 or 2", s
IF STYP(s) == SCOD("BINARY_1")
   GETEXTRADATA 1, s
   maxterm = VEC1(1)
   ord = 0
   sum = 0
   LABEL 10
      ord = ord + 1
      sum = sum + ord + 1
   IF sum < maxterm THEN GOTO 10
   LABEL 11
   IF maxterm < sum
      maxterm=maxterm+1
      VEC1(maxterm+2) = 0.0
      GOTO 11
   ENDIF  
   r0 = VEC1(2)
   INPUT "Enter x coordinate of y-scan", x0
   INPUT "Enter y coordinate of x-scan", y0
   PRINT "diffractive zone width for surface  ",s
   PRINT "y direction scan at x=",x0
   PRINT "y coord ",un$,"        zone width (microns)"
   x = x0/r0
   FOR y1=SDIA(s),-SDIA(s),-SDIA(s)/10
      y = y1/r0
      GOSUB XDeriv
      GOSUB YDeriv
      dp = SQRT(dpdx*dpdx + dpdy*dpdy)
      IF dp>1.0E-6
         wid = const*r0/dp
         PRINT y1,"                ", wid
      ELSE
         PRINT y1,"                    inf"
      ENDIF
   NEXT
   PRINT "x direction scan at y=",y0
   PRINT "x coord ",un$,"        zone width (microns)"
   y = y0/r0
   FOR x1=SDIA(s),-SDIA(s),-SDIA(s)/10
      x = x1/r0
      GOSUB XDeriv
      GOSUB YDeriv
      dp = SQRT(dpdx*dpdx + dpdy*dpdy)
      IF dp>1.0E-6
         wid = const*r0/dp
         PRINT x1,"                ", wid
      ELSE
         PRINT x1,"                    inf"
      ENDIF
   NEXT
ENDIF
IF STYP(S) == SCOD("BINARY_2")
   GETEXTRADATA 1, s
   r0 = VEC1(2)
   PRINT "diffractive zone width for surface  ",s
   PRINT "radius ",un$,"        zone width (microns)"
   FOR r1=SDIA(s),0,-SDIA(s)/10
      r = r1/r0
      GOSUB RDeriv
      dp = ABSO(dpdr)
      IF dp>1.0E-6
         wid = const*r0/dp
         PRINT r1,"                ", wid
      ELSE
         PRINT r1,"                    inf"
      ENDIF
   NEXT
ENDIF
END
!------------------------------------------------
! Subroutine for partial derivative of phi wrt x
SUB XDeriv
   dphi = 0.0
   i = ord
   LABEL 1
   IF i>1
      row = 0.0
      j = ord - i
      LABEL 2
      IF j>0
         k=(i+j)*(i+j+1)/2 + j
         row = (row + VEC1(k+2))*y
         j=j-1
         GOTO 2
      ENDIF
      k = i*(i+1)/2
      row = (row + VEC1(k+2))*i
      dphi = (dphi + row)*x
      i=i-1
      GOTO 1      
   ENDIF
   row = 0.0
   j=ord-1
   LABEL 3
   IF j>0
      k = (j+1)*(j+2)/2 + j
      row =(row + VEC1(k+2))*y
      j=j-1
      GOTO 3
   ENDIF
   dpdx = dphi + row + VEC1(3)
RETURN
!------------------------------------------------
! Subroutine for partial derivative of phi wrt y
SUB YDeriv
   dphi = 0.0
   j=ord
   LABEL 4
   IF j>1
      col = 0.0
      i=ord-j
      LABEL 5
      IF i>0
         k = (i+j)*(i+j+1)/2 + j
         col = (col + VEC1(k+2))*x
         i=i-1
         GOTO 5
      ENDIF
      k = j*(j+1)/2 + j
      col = (col + VEC1(k+2))*j
      dphi = (dphi + col)*y
      j=j-1
      GOTO 4
   ENDIF
   col = 0.0
   i=ord-1
   LABEL 6
   IF i>0
      k = (i+1)*(i+2)/2 + 1
      col = (col + VEC1(k+2))*x
      i=i-1
      GOTO 6
   ENDIF
   dpdy = dphi + col + VEC1(4)
RETURN
!-------------------------------------------
! Subroutine for derivative of phi wrt to r
SUB RDeriv
   dpdr = 0.0
   r2 = r*r
   FOR j=VEC1(1),1,-1
      dpdr = dpdr*r2 + 2*j*VEC1(j+2)*r 
   NEXT
RETURN
!--------------------------------------------
! Subroutine to test above code, not used otherwise
SUB Test
Input "Enter hx", hx
Input "Enter hy", hy
Input "Enter px", px
Input "Enter py", py
PRINT "diffractive zone width for surface  ",s,"  based on ray trace"
RAYTRACE hx,hy,px,py,PWAV()
dx = ABSO( INDX(s) * RAYL(s) - INDX(s-1) * RAYL(s-1) )
dy = ABSO( INDX(s) * RAYM(s) - INDX(s-1) * RAYM(s-1) )
d= SQRT( dx*dx + dy*dy )
IF d != 0
   PRINT "X         Y             zone size"
   PRINT RAYX(s),"   ",RAYY(s),"        ", WAVL( PWAV() ) / d
ENDIF
x1 = RAYX(s)
y1 = RAYY(s)
r1 = SQRT(x1*x1 + y1*y1)
RETURN