'---------------------------------------------------------------
'                 DIFFRACT.BAS
'
'  http://www.achilles.net/~jtalbot/amateur/diffract
'
'  by John Talbot, jtalbot@achilles.net  Apr 19,1997.
'
'  This version runs only on a PC 
'
' Computes the fresnel diffration pattern of various obstacles
' and plots them to the screen as graphs and/or intensity coded
' colors. The pattern can also be plotted to an HPGL/2 
' LaserJetIII printer file.
'
'(Minor Bug : The gaussian waist laser cross section (waist)
' isn't re-computed when the detector distance changes.
' Should modify program to change the quadratic exponential
' to take this into account.)
'
'----------------------------------------------------------------
 DEFDBL A-Z
 DECLARE SUB Fresnel (C#, S#, XOrigin#, XPos#)
 DECLARE SUB Intensity (X#, Y#, Z#)
 DECLARE SUB Fresnel3 (C#, S#, X#)
 DECLARE SUB MapXY (X#, Y#, X$, Y$)
 DECLARE SUB VideoXY (X#, Y#, XX#, YY#)
 DECLARE SUB LaserXY (X#, Y#, X$, Y$)
 DECLARE SUB Vector (X1#, Y1#, X2#, Y2#, TheColor#)
 DECLARE SUB LaserVector (X1#, Y1#, X2#, Y2#, ThePenWidth#)
 DECLARE SUB GetNumber (Number#)
 DECLARE SUB GetKey (Choice#, Maximum#)
 DECLARE SUB BesselJ1 (Rho#, j1#)

 DIM SHARED Y!(3300), PI, A, B, j(500), B1(11), B2(11), B3(11), B4(11)
 DIM SHARED XBorderL, YBorderL, XScaleL, YScaleL, DomainMin, DomainMax, Min, Max
 DIM SHARED XBorder, YBorder, XScale, YScale, ScreenY, Plotter, MaxPt, XRange

 GOSUB Constants
 GOSUB Default

 WHILE Plot <> Quit
   GOSUB MenuChoice
   X = DomainMin
   dX = XRange / (MaxPt - 1)
   FOR Pt = 0 TO MaxPt - 1
     GOSUB Amplitude
     GOSUB Convolve
     Y!(Pt) = YJ     'Function of X
     X = X + dX
   NEXT Pt

   IF Plot = Laser OR Plot = Video THEN
     GOSUB LaserPlot
   ELSEIF Plot = TellAGraf THEN
     GOSUB TellAGrafFile
   ELSEIF Plot = Simulation THEN
     IF Obstacle = Rectangle THEN
       GOSUB Rectangular
     ELSE
       GOSUB LaserSpot
     END IF
   END IF
   WHILE INKEY$ = "": WEND
 WEND

END

'---------------------------------------------------------------------
Default:

  'default parameters   all units are in mm (millimeters)
  DomainMin = -20  'Start scan at detector  (usually -80 or -30)
  DomainMax = 20   'Stop scan at detector   (usually +80 or +30)

  WG = 52.1        'Gaussian envellope width at 13.1% (for fixed A1+A2)
  A1 = 88          'Source to Obstacle distance
  A2 = 1469        'Obstacle to Screen distance
  Wavelength = 6.3282E-04  'Laser wavelength;  632.82 nm (monomode only)
  Diameter = .13   'IF Obstacle = Wire      then Diameter
                   'IF Obstacle = Slit      then Width
                   'IF Obstacle = Rectangle then Width
  Height = .25     'IF Obstacle = Rectangle then Aperture Height
  FirstSlit = .1   'IF Obstacle = TwoSlit   then First slit width
  SecondSlit = .1  'IF Obstacle = TwoSlit   then Second slit width
  SlitSpacing = .2 'IF Obstacle = TwoSlit   then Slit separation
  EdgeCenter = 0   'IF Obstacle = Edge      then Center of Edge

  GaussWidth = -WG * WG / 2   'should compute this value
  BesselZero = .8 * WG        'rough initial estimate of first zero
  OverlayFlag = 1             'if 1 then plot enveloppe overlay

RETURN

'------------------------------------------------------------------------
MenuChoice:

 EdgeColor = 10
 TitleColor = 14
 BorderColor = 5
 OverlayColor = 6
 PlotColor = 15
 SCREEN 12: CLS : COLOR PlotColor
 PRINT "     Fresnel Diffraction          by John Talbot, Feb 8,1993"
 PRINT "     ------------------- "
 PRINT
 PRINT " Computer simulation of diffraction patterns produced by various "
 PRINT " simple obstacles using the Fresnel approximation."
 PRINT " (ref: Miles V. Klein 'Optics',1970 p.363-392, and AJP 40(1972)74)"
 PRINT
 PRINT "    1. None"
 PRINT "    2. WIRE"
 PRINT "    3. Edge"
 PRINT "    4. Slit"
 PRINT "    5. Two slits"
 PRINT "    6. Rectangular aperture"
 PRINT "    7. Quit"
 PRINT " What type of obstacle ?";
 GetKey Obstacle, 7#
 IF Obstacle = Quit THEN END
 PRINT
 PRINT "    1. Uniform     Isotropic Point Source (theoretical)"
 PRINT "    2. GAUSSIAN    Spatially filtered expanding beam"
 PRINT "    3. Bessel      Laser improperly focussed on pinhole"
 PRINT " What type of laser beam (angular envellope) ?";
 GetKey Beam, 3#
 PRINT
 IF Obstacle <> Rectangle THEN
  PRINT "    1. Screen graph of intensity versus position"
  PRINT "    2. Simulated laser spot pattern (in 2 dimensions)"
  PRINT "    3. Animation sweeping one parameter in time"
  PRINT "    4. Laser printer    (LaserJetIIIp HPGL/2)"
  PRINT "    5. TellAGraf script"
  PRINT " What type of plot ?";
  GetKey Plot, 5#
  PRINT
  IF Plot = Laser OR Plot = TellAGraf THEN
    Plotter = Plotter + 1
    Plotter$ = STR$(Plotter): FileName$ = CurrentDir$ + "Plot" + RIGHT$(Plotter$, LEN(Plotter$) - 1)
    PRINT "  Will be saved to disk as "; FileName$
    PRINT
  END IF
  IF Plot = Laser OR Plot = Video THEN
     PRINT "    1. Plot overlay of undisturbed laser intensity"
     PRINT "    2. No Overlay"
     PRINT " Plot Characteristics ?";
     GetKey OverlayFlag, 2#
     PRINT
  END IF

 ELSE
  Plot = Simulation
 END IF
 PRINT "  Press Return if you don't want to change the previous value ..."
 PRINT USING "  #.####^^^^ mm  What laser wavelength"; Wavelength; : GetNumber Wavelength
 PRINT USING "  ####.# mm      Source to obstacle distance"; A1; : GetNumber A1
 PRINT USING "  ####.# mm      Obstacle to screen distance"; A2; : GetNumber A2
 A = SQR(2 * A1 / (Wavelength * A2 * (A1 + A2)))  'Scaling down and SQR(2/(Wavelength*L))
 B = (A1 + A2) / A1                       'Scaling up
 
 PRINT USING "  ####.# mm      Detector starts scanning at"; DomainMin; : GetNumber DomainMin
 PRINT USING "  ####.# mm      Detector stops scanning at"; DomainMax; : GetNumber DomainMax
 XRange = DomainMax - DomainMin
 IF XRange <= 0 THEN PRINT "Illegal Range ": END

 IF Obstacle = Wire THEN
   PRINT USING "  ##.### mm      Wire diameter"; Diameter; : GetNumber Diameter
   Title$ = "FRESNEL DIFFRACTION PATTERN OF A WIRE"
   XB = Diameter / 2: XA = -XB
 ELSEIF Obstacle = Slit THEN
   PRINT USING "  ##.### mm      Slit width"; Diameter; : GetNumber Diameter
   Title$ = "FRESNEL DIFFRACTION PATTERN OF A SLIT"
   XB = Diameter / 2: XA = -XB
 ELSEIF Obstacle = TwoSlit THEN
   PRINT USING "  ##.### mm      First slit width"; FirstSlit; : GetNumber FirstSlit
   PRINT USING "  ##.### mm      Second slit width"; SecondSlit; : GetNumber SecondSlit
   PRINT USING "  ##.### mm      Slit separation"; SlitSpacing; : GetNumber SlitSpacing
   Title$ = "FRESNEL DIFFRACTION PATTERN OF TWO SLITS"
   XX = SlitSpacing / 2
   XZ = -SlitSpacing / 2
   XB = SlitSpacing / 2 + SecondSlit
   XA = -SlitSpacing / 2 - FirstSlit
 ELSEIF Obstacle = Edge THEN
   PRINT USING "  ##.### mm      Obstacle Center"; EdgeCenter; : GetNumber EdgeCenter
   Title$ = "FRESNEL DIFFRACTION PATTERN OF A STRAIGHT EDGE"
   XA = EdgeCenter
 ELSEIF Obstacle = Rectangle THEN
   PRINT USING "  ##.### mm      Aperture width"; Diameter; : GetNumber Diameter
   PRINT USING "  ##.### mm      Aperture height"; Height; : GetNumber Height
   Title$ = "FRESNEL DIFFRACTION PATTERN OF RECTANGULAR APERTURE"
   XB = Diameter / 2: XA = -XB
 ELSEIF Obstacle = None THEN
   Title$ = "UNDISTURBED ANGULAR PATTERN OF LASER (NO OBSTACLE)"
 END IF

 IF Plot = Simulation THEN
   MAXX = 320: MAXY = 200
   MaxPt = MAXX
   Aspect = 1.3
   HeadRoom = 0    ' 0% beyond extremums
   XBorder = 0
   YBorder = 0
   ScreenX = MAXX - 2 * XBorder
   ScreenY = MAXY - 2 * YBorder
   EdgeColor = 65
   TitleColor = 64
   BorderColor = 66
   OverlayColor = 64
   PlotColor = 63
   GOSUB GreyScales
 ELSE
   MAXX = 640: MAXY = 480
   HeadRoom = 2 / 100        ' 2% beyond extremums
   XBorder = 5
   YBorder = 18
   ScreenX = MAXX - 2 * XBorder
   ScreenY = MAXY - 2 * YBorder
 END IF

 IF Plot = Video OR Plot = Animation THEN
   MaxPt = ScreenX
 ELSEIF Plot = Laser THEN
   MaxPt = ScreenXL  'DO NOT not change this!
 ELSEIF Plot = TellAGraf THEN
   MaxPt = ScreenXL  'can be changed
 END IF

RETURN

'-----------------------------------------------------------------
Constants:

  DATA  1.595769140, -0.000000033,  0.000000000,  0.199471140
  DATA -0.000001702,  4.255387524, -0.024933975,  0.000000023
  DATA -6.808568854, -0.000092810,  0.000003936, -0.009351341
  DATA -0.000576361, -7.780020400,  0.005770956,  0.000023006
  DATA  6.920691902, -0.009520895,  0.000689892,  0.004851466
  DATA -0.016898657,  5.075161298, -0.009497136,  0.001903218
  DATA -3.050485660, -0.138341947,  0.011948809, -0.017122914
  DATA -0.075752419, -1.363729124, -0.006748873,  0.029064067
  DATA  0.850663781, -0.403349276,  0.000246420, -0.027928955
  DATA -0.025639041,  0.702222016,  0.002102967,  0.016497308
  DATA -0.150230960, -0.216195929, -0.001217930, -0.005598515
  DATA  0.034404779,  0.019547031,  0.000233939,  0.000838386

  RESTORE Constants
  FOR i = 0 TO 11
    READ B1(i), B2(i), B3(i), B4(i)
  NEXT i
  
  None = 1
  Wire = 2
  Edge = 3
  Slit = 4
  TwoSlit = 5
  Rectangle = 6
  Quit = 7: Obstacle = Wire          'default type of obstacle

  Video = 1
  Simulation = 2
  Animation = 3
  Laser = 4
  TellAGraf = 5: Plot = Video       'default type of plot

  Uniform = 1
  Gaussian = 2
  Bessel = 3: Beam = Gaussian        'default type beam

  PI = 3.141592653589#

  'There must be a bug here! (improper scaling somewhere!)

  MAXX = 640: MAXY = 480       'Default Screen resolution
  MAXXL = 3300: MAXYL = 2400   'Laser printer resolution (HP LaserJetIIIp)

  XBorderL = MAXXL * 8 / 100   'should leave room for axis labels
  YBorderL = MAXYL * 8 / 100   'and for axis coordinates
  ScreenXL = MAXXL - 2 * XBorderL
  ScreenYL = MAXYL - 2 * YBorderL
  'Make shure that MaxPt=ScreenXL otherwise PR (HPGL/2) command won't work

  XAxis$ = "POSITION (mm)"
  YAxis$ = "INTENSITY (arbitrary units)"

  PlotPen = .09     'mm or  2 dots width
  TextPen = .17     'mm or  2 dots width
  TitlePen = .34    'mm or  4 dots width
  BorderPen = .34   'mm or  4 dots width
  ThickPen = .51    'mm or  6 dots width
  FatPen = 1.19     'mm or 14 dots width
  mm = 25.4 / 300   'Multiply Dots by mm to obtain millimeters (300 DPI)
  PenSize = TextPen

  CurrentDir$ = "C:"        'Output directory
  Plotter = 0        'PCL5 and HPGL/2 plot file number for LaserJetIIIp.

RETURN

'----------------------------------------------------------------
'Set the correct integration limits and
'compute the Intesity of the diffraction patern.

Amplitude:
    IF Obstacle = Wire THEN
      Fresnel C1, S1, XA, X
      Fresnel C2, S2, XB, X
      Intensity 1# - C2 + C1, 1# - S2 + S1, YJ
    ELSEIF Obstacle = Edge THEN
      Fresnel C1, S1, XA, X
      Intensity .5# - C1, .5# - S1, YJ
    ELSEIF Obstacle = Slit OR Obstacle = Rectangle THEN
      Fresnel C1, S1, XA, X
      Fresnel C2, S2, XB, X
      Intensity C2 - C1, S2 - S1, YJ
    ELSEIF Obstacle = TwoSlit THEN
      Fresnel C1, S1, XA, X
      Fresnel C2, S2, XB, X
      Fresnel C3, S3, XX, X
      Fresnel C4, S4, XZ, X
      Intensity C2 - C1 - C3 + C4, S2 - S1 - S3 + S4, YJ
    ELSEIF Obstacle = None THEN
      YJ = 1
    END IF
RETURN
'-------------------------------------------------------------------
'Multiply the diffraction pattern by the Laser Beam Angular Envellope

Convolve:
    IF Beam = Gaussian THEN
       YJ = 10 * YJ * EXP(X * X / GaussWidth)
    ELSEIF Beam = Bessel AND Plot <> Simulation THEN
       Rho = 3.832 * X / BesselZero
       BesselJ1 Rho, j1
       j1 = 2 * j1 / Rho
       YJ = 10 * YJ * j1 * j1
    ELSEIF Beam = Uniform THEN
       YJ = 10 * YJ
    END IF
RETURN

'-------------------------------------------------------------------------
'Laser Printer plotting and screen previewer
' (after leaving Basic type:  COPY plot1 LPT1: , COPY plot2 LPT1:  etc...)

LaserPlot:
   CLS
   GOSUB Scale
   GOSUB Coordinates
   GOSUB Shadow
   IF Plot = Laser THEN
     GOSUB InitPrinter
     GOSUB LaserCoordinates
     GOSUB LaserShadow
   END IF
   IF OverlayFlag = 1 AND Obstacle <> None THEN GOSUB OverLay

   COLOR TitleColor: LOCATE 1, 5: PRINT Title$
   LOCATE 30, 1: PRINT USING "###.##"; DomainMin;
   LOCATE 30, 74: PRINT USING "###.##"; DomainMax;
   LOCATE 30, 35: PRINT XAxis$; : LOCATE 1, 1:
   X = DomainMin
   dX = XRange / (MaxPt - 1)
   Y = Y!(0)
   VideoXY X, Y, XS, YS
   PSET (XS, YS), PlotColor

   'Start first point on plot
   Yp = INT(YBorderL + YScaleL * (Y!(0) - Min))
   IF Plot = Laser THEN
     Y = Y!(0): LaserXY X, Y, X$, Y$
     PRINT #Plotter, USING "PU; PW ##.##;"; PlotPen;
     PRINT #Plotter, "PA"; X$; Y$; "; PD;"
   END IF
   Y$ = "PR"
   FOR j = 1 TO MaxPt - 1
      X = X + dX
      IF Plot = Laser THEN
        Y = INT(YBorderL + YScaleL * (Y!(j) - Min))
        Y$ = Y$ + " 1" + STR$(Y - Yp)
        Yp = Y
        IF LEN(Y$) > 75 THEN PRINT #Plotter, Y$; ";": Y$ = "PR"
      END IF
      Y = Y!(j): VideoXY X, Y, XS, YS
      LINE -(XS, YS), PlotColor
   NEXT j
   IF LEN(Y$) > 2 AND Plot = Laser THEN PRINT #Plotter, Y$; ";"

   IF Plot = Laser THEN
     'Restore defaults, Enter PCL Mode, Force page out and reset the printer
     PRINT #Plotter, "SP1; LT; TR1;"; CHR$(27); "%1A"; CHR$(27); "E"
     CLOSE Plotter
   END IF

RETURN
'----------------------------------------------------------------------
'Plot on top of diffraction pattern the undisturbed laser beam

OverLay:

   X = DomainMin
   dX = XRange / (MaxPt - 1)
   YJ = 1: GOSUB Convolve  'Output is YJ
   VideoXY X, YJ, XS, YS
   PSET (XS, YS), OverlayColor

   'Start first point on plot of undisturbed beam
   Yp = INT(YBorderL + YScaleL * (YJ - Min))
   IF Plot = Laser THEN
     LaserXY X, YJ, X$, Y$
     PRINT #Plotter, USING "PU; LT2,2,0; PW ##.##;"; PlotPen;
     PRINT #Plotter, "PA"; X$; Y$; "; PD;"
   END IF
   Y$ = "PR"
   FOR j = 1 TO MaxPt - 1
      X = X + dX
      YJ = 1: GOSUB Convolve
      IF Plot = Laser THEN
        Y = INT(YBorderL + YScaleL * (YJ - Min))
        Y$ = Y$ + " 1" + STR$(Y - Yp)
        Yp = Y
        IF LEN(Y$) > 75 THEN PRINT #Plotter, Y$; ";": Y$ = "PR"
      END IF
      VideoXY X, YJ, XS, YS
      LINE -(XS, YS), OverlayColor
   NEXT j
   IF LEN(Y$) > 2 AND Plot = Laser THEN PRINT #Plotter, Y$; ";"
   IF Plot = Laser THEN PRINT #Plotter, "LT;"
RETURN

'------------------------------------------------------
LaserCoordinates:

   'Print axis labels and title

   PRINT #Plotter, USING "PU; PW ##.##;"; TextPen;
   PRINT #Plotter, "DT#; SB; SD; SR 1.1,2.0; LO 14; DR 0,1; PA";
   PRINT #Plotter, STR$(INT(XBorderL / 2)); STR$(INT(MAXYL / 2)); ";PD;";
   PRINT #Plotter, "LB"; YAxis$; "#; PU;"

   PRINT #Plotter, "LO 16; DR 1,0; PA";
   PRINT #Plotter, STR$(INT(MAXXL / 2)); STR$(INT(YBorderL / 2)); "; PD;";
   PRINT #Plotter, "LB"; XAxis$; "#; PU;"

   PRINT #Plotter, USING "PW ##.##;"; TitlePen;
   PRINT #Plotter, "LO 14; DR 1,0; PA";
   PRINT #Plotter, STR$(INT(MAXXL / 2)); STR$(INT(MAXYL - YBorderL)); "; PD;";
   PRINT #Plotter, "LB"; Title$; "#; PU;"

   'Draw Box around graph

   LaserVector DomainMin, Min, DomainMax, Min, BorderPen
   LaserVector DomainMax, Min, DomainMax, Max, BorderPen
   LaserVector DomainMax, Max, DomainMin, Max, BorderPen
   LaserVector DomainMin, Max, DomainMin, Min, BorderPen

   'Print axis ticks and values
   'XAxis
   XTicSize = 10: IF XRange > 100 THEN XTicSize = XTicSize * 2
   k = INT(DomainMin / XTicSize) + 1
   PRINT #Plotter, "LO16; DR 1,0; "
   WHILE k * XTicSize < DomainMax
     X = k * XTicSize
     X$ = STR$(X): IF X >= 0 THEN X$ = RIGHT$(X$, LEN(X$) - 1)
     LaserVector X, Min, X, Min + YRange * .015#, BorderPen
     LaserXY X, Min, X1$, Y1$
     PRINT #Plotter, USING "PW ##.##; "; TextPen;
     PRINT #Plotter, "PA"; X1$; Y1$; "; LB"; X$; "#; "
     k = k + 1
   WEND

   YTicSize = 2
   k = INT(Min / YTicSize) + 1
   PRINT #Plotter, "LO18; "
   WHILE k * YTicSize < Max
     Y = k * YTicSize
     Y$ = STR$(Y): IF Y >= 0 THEN Y$ = RIGHT$(Y$, LEN(Y$) - 1)
     LaserVector DomainMin, Y, DomainMin + XRange * .015#, Y, BorderPen
     LaserXY DomainMin, Y, X1$, Y1$
     PRINT #Plotter, USING "PW ##.##; "; TextPen;
     PRINT #Plotter, "PA"; X1$; Y1$; "; LB"; Y$; "#; "
     k = k + 1
   WEND
RETURN
'----------------------------------------------------------------------
LaserShadow:
  Y = Min + YRange * .03
  'LaserVector automatically performs hardclip if out of bounds
  IF Obstacle = Wire OR Obstacle = Slit THEN
    LaserVector B * XA, Y, B * XB, Y, FatPen
  ELSEIF Obstacle = Edge THEN
    LaserVector DomainMin, Y, B * XA, Y, FatPen
  ELSEIF Obstacle = TwoSlit THEN
    LaserVector B * XA, Y, B * XZ, Y, FatPen
    LaserVector B * XX, Y, B * XB, Y, FatPen
  END IF
RETURN
'-----------------------------------------------------
'Plot Geometric shadow Edges of Square aperture
Shadow:
  IF Obstacle = Wire OR Obstacle = Slit THEN
    Vector B * XA, Min, B * XA, Max, EdgeColor
    Vector B * XB, Min, B * XB, Max, EdgeColor
  ELSEIF Obstacle = Edge THEN
    Vector B * XA, Min, B * XA, Max, EdgeColor
  ELSEIF Obstacle = TwoSlit THEN
    Vector B * XA, Min, B * XA, Max, EdgeColor
    Vector B * XX, Min, B * XX, Max, EdgeColor
    Vector B * XZ, Min, B * XZ, Max, EdgeColor
    Vector B * XB, Min, B * XB, Max, EdgeColor
  END IF
RETURN

Coordinates:
   'Draw Box around graph
   VideoXY DomainMin, Min, X1, Y1
   VideoXY DomainMax, Max, X2, Y2
   LINE (X1, Y1)-(X2, Y2), BorderColor, B
RETURN
'------------------------------------------------------
InitPrinter:

 OPEN FileName$ FOR OUTPUT AS Plotter

 '----LaserJet IIIp Plotter commands---
 ' Reset Printer
 ' PageSize=Letter, Orientation=Portrait, TopMargin=0
 ' PCL CursorX=0 dots, PCL CursorY=0 dots
 ' Set picture frame to 8 X 11 inch (decipoints), Anchor point at PCL cursor
 ' Enter HPGL/2
 ' Initialize HPGL/2 environment (Origin at bottom left)
 ' (5760 decipoints = 8 inch, 7920 decipoints = 11 inch )
 ' (2400 pixels     = 8 inch, 3300 pixels     = 11 inch )
 ' ( 720 decipoints = 1 inch,  2.4 decipoints =  1 pixel)

 PRINT #Plotter, CHR$(27); "E";
 PRINT #Plotter, CHR$(27); "&l2a1o0E";
 PRINT #Plotter, CHR$(27); "*p0x0Y";
 PRINT #Plotter, CHR$(27); "*c7920x5760y0T";
 PRINT #Plotter, CHR$(27); "%1B"
 PRINT #Plotter, "IN;"

 ' Select Pen 1
 ' Input Relative control points at bottom left and top right(normal),
 ' Set Custom user coordinates to 1/300 inch units
 ' ( 2400 dot units = 8 inch , 3300 dot units = 11 inch )
 ' Line Atributes : line ends = round, line joins = round
 ' default Line Type
 PRINT #Plotter, "SP1; SC 0,3300,0,2400; LA 1,4,2,4; LT;"

RETURN

'--------------------------------------------
' all variable ending in 'L' are for the Laser Printer

Scale:
    Max = -1E+20: Min = 1E+20
    FOR Pt = 0 TO MaxPt - 1
      Y = Y!(Pt)
      IF Y > Max THEN Max = Y
      IF Y < Min THEN Min = Y
    NEXT Pt

    YRange = Max - Min
    IF YRange = 0 THEN YRange = Min / 10
    IF YRange = 0 THEN YRange = 1
    Min = Min - YRange * HeadRoom
    Max = Max + YRange * HeadRoom

    YRange = Max - Min
    IF YRange = 0 THEN YRange = 1

    YScale = ScreenY / YRange: YScaleL = ScreenYL / YRange

    XRange = DomainMax - DomainMin
    XScale = ScreenX / XRange: XScaleL = ScreenXL / XRange
RETURN
'----------------------------------------------------------

TellAGrafFile:
  OPEN FileName$ FOR OUTPUT AS Plotter
  q$ = CHR$(39)
  PRINT #Plotter, "GENERATE A PLOT."
  PRINT #Plotter, "X AXIS LABEL TEXT IS " + q$ + XAxis$ + q$ + "."
  PRINT #Plotter, "Y AXIS LABEL TEXT IS " + q$ + YAxis$ + q$ + "."
  PRINT #Plotter, "TITLE TEXT IS " + q$ + Title$ + q$ + "."
  PRINT #Plotter, "INPUT DATA."

  X = DomainMin
  dX = XRange / (MaxPt - 1)
  FOR Pt = 0 TO MaxPt - 1
    PRINT #Plotter, USING "####.#### ##.####"; X; Y!(Pt)
    X = X + dX
  NEXT Pt
  PRINT #Plotter, "END OF DATA."
  PRINT #Plotter, "GO."
  CLOSE Plotter
RETURN



'------------------------------------------------------------
'------------------------------------------------------------
'Color Graphics Intensive Section of program
'------------------------------------------------------------
GreyScales:

   'Approximates laser intensity using 64 levels of quantization,
   ' and roughly simulates laser color according to Wavelength.
   SCREEN 13
   IF Wavelength * 1000000# > 620 THEN 'RED
     Red = 1: Green = 0: Blue = 0
   ELSEIF Wavelength * 1000000# > 570 THEN 'YELLOW
     Red = 1: Green = 256&: Blue = 0
   ELSEIF Wavelength * 1000000# > 520 THEN 'GREEN
     Red = 0: Green = 256&: Blue = 0
   ELSEIF Wavelength * 1000000# > 480 THEN 'CYAN
     Red = 0: Green = 256&: Blue = 256& * 256&
   ELSE 'BLUE
   Red = 0: Green = 0: Blue = 256& * 256&
   END IF

   'This method minimizes color fringing by linearly increasing
   'the intensity by integer steps for each color component separately)
   FOR i = 0 TO 63
      RGBColor& = Red * i + Green * i + Blue * i
      PALETTE i, RGBColor&
      LINE (i * ScreenX / 64, ScreenY / 3)-STEP(ScreenX / 64 - 1, ScreenY / 3), i, BF
   NEXT i
   RGBColor& = 63: PALETTE 64, RGBColor&
   RGBColor& = 256 * 63: PALETTE 65, RGBColor&
   RGBColor& = 256& * 256& * 63&: PALETTE 66, RGBColor&
   RGBColor& = 63& + 256& * 63&: PALETTE 67, RGBColor&

   LINE (0, ScreenY / 3)-(ScreenX - 1, 2 * ScreenY / 3), 63, B
   COLOR 63: LOCATE 6, 1: PRINT "Plotting in Progress ...."
   COLOR 66: LOCATE 20, 5: PRINT "Adjust your monitor for linear scale"
   COLOR 65: LOCATE 22, 1: PRINT Title$; : LOCATE 1, 1

RETURN

'------------------------------------------------------------
'Plot the diffraction in 2 dimensions of a rectangular hole

Rectangular:
   GOSUB Scale
   MaxSpot = Max - YRange * HeadRoom   'reduce this for fainter detail

   Min = -Aspect * ScreenY / (2 * XScale)
   Max = -Min
   YRange = Max - Min
   YScale = ScreenY / YRange
   X1 = B * XA: X2 = B * XB
   XA = Height / 2: XB = -XA
   Y1 = B * XA: Y2 = B * XB
  
   FOR j = 0 TO ScreenY - 1
    Y = Aspect * (j - ScreenY / 2) / XScale     'in millimeters
    X = Y: GOSUB Amplitude
    YJ = 64 * YJ * EXP(Y * Y / GaussWidth) / MaxSpot
    FOR i = 0 TO ScreenX - 1
     X = INT(i * (MaxPt - 1) / (ScreenX - 1))
     Z = YJ * Y!(X)
     IF Z > 63 THEN Z = 63     'set to maximum saturation
     PSET (i, j), Z
    NEXT i
   NEXT j
   'LOCATE 23, 1: PRINT "Type any key to Return to Menu"; : LOCATE 1, 1
   WHILE INKEY$ = "": WEND
   GOSUB ApertureEdge
RETURN

ApertureEdge:
   'Plot Geometric shadow Edges of Square aperture
   Vector X1, Y1, X2, Y1, EdgeColor
   Vector X2, Y1, X2, Y2, EdgeColor
   Vector X2, Y2, X1, Y2, EdgeColor
   Vector X1, Y2, X1, Y1, EdgeColor
RETURN

'------------------------------------------------------
LaserSpot:
   GOSUB Scale
   MaxSpot = Max - YRange * HeadRoom   'reduce this for fainter detail
   IF Beam = Gaussian THEN
    FOR j = 0 TO ScreenY - 1
     Y = Aspect * (j - MAXY / 2) / XScale       'in millimeters
     YJ = 64 * EXP(Y * Y / GaussWidth) / MaxSpot
     FOR i = 0 TO ScreenX - 1
      X = INT(i * (MaxPt - 1) / (MAXX - 1))
      Z = YJ * Y!(X)
      IF Z > 63 THEN Z = 63     'set to maximum saturation
      PSET (i, j), Z
     NEXT i
    NEXT j
   ELSEIF Beam = Uniform THEN
    FOR i = 0 TO ScreenX - 1
      X = INT(i * (MaxPt - 1) / (ScreenX - 1))
      Z = 64 * Y!(X) / MaxSpot
      IF Z > 63 THEN Z = 63     'set to maximum saturation
      LINE (i, 0)-(i, MAXY - 1), Z
     NEXT i
   ELSEIF Beam = Bessel THEN
    FOR j = 0 TO ScreenY - 1
     Y = Aspect * (j - ScreenY / 2) / XScale         'in millimeters
     FOR i = 0 TO ScreenX - 1
      X = INT(i * (MaxPt - 1) / (ScreenX - 1))
      XX = DomainMin + i * XRange / (ScreenX - 1)    'in millimeters
      Z = Y!(X) / MaxSpot
      Rho = 3.832 * SQR(XX * XX + Y * Y) / BesselZero
      BesselJ1 Rho, j1
      j1 = 2 * j1 / Rho
      Z = 64 * Z * j1 * j1
      IF Z > 63 THEN Z = 63     'set to maximum saturation
      PSET (i, j), Z
     NEXT i
    NEXT j
   END IF

   WHILE INKEY$ = "": WEND
   X = DomainMin: dX = XRange / (MaxPt - 1)
   Temp = Y!(0)
   VideoXY X, Temp, X1, Y1
   PSET (X1, Y1), PlotColor
   FOR i = 0 TO ScreenX - 1
      Temp = Y!(i)
      VideoXY X, Temp, X1, Y1
      LINE -(X1, Y1), PlotColor
      X = X + dX
   NEXT i
   GOSUB Shadow

   
RETURN

PlotEdge:

RETURN

'----------------------------------------------------------
'Numerical evaluation of Fresnel integrals using the Tau-method
'of Lanczos. J.Boersma, 'Mathematics of Computation' 14(1960)380
'Accurate to better than 2 X 10^-9

SUB BesselJ1 (Rho, j1) STATIC
 n% = 1.1 * ABS(Rho) + 15
 IF n% > 499 THEN n% = 499
 j(n%) = 0
 j(n% - 1) = 1E-12
 Factor = 2 / Rho
 FOR k% = n% - 1 TO 1 STEP -1
   j(k% - 1) = (Factor * k%) * j(k%) - j(k% + 1)
 NEXT k%
 Sum = 0
 FOR k% = 2 TO n% STEP 2
   Sum = Sum + j(k%)
 NEXT k%
 Sum = j(0) + 2 * Sum
 j1 = j(1) / Sum
END SUB

SUB Fresnel (C, S, XOrigin, XPos) STATIC
  X = A * (B * XOrigin - XPos)
  Fresnel3 C, S, X
END SUB

SUB Fresnel3 (C, S, X) STATIC
  U = ABS(X * X * PI / 2#)
  IF U > 4 THEN
    Z = 4 / U
    Sum3 = 0
    Sum4 = 0
    FOR i% = 11 TO 1 STEP -1
      Sum3 = (Sum3 + B3(i%)) * Z
      Sum4 = (Sum4 + B4(i%)) * Z
    NEXT i%
    Sum3 = Sum3 + B3(0)
    Sum4 = Sum4 + B4(0)
    CosX = COS(U): SinX = SIN(U)
    Z = SQR(Z)
    C = .5 + Z * (CosX * Sum3 + SinX * Sum4)
    S = .5 - Z * (CosX * Sum4 - SinX * Sum3)
  ELSEIF U > 0 THEN
    Z = U / 4
    Sum1 = 0
    Sum2 = 0
    FOR i% = 11 TO 1 STEP -1
      Sum1 = (Sum1 + B1(i%)) * Z
      Sum2 = (Sum2 + B2(i%)) * Z
    NEXT i%
    Sum1 = Sum1 + B1(0)
    Sum2 = Sum2 + B2(0)
    CosX = COS(U): SinX = SIN(U)
    Z = SQR(Z)
    C = Z * (CosX * Sum1 + SinX * Sum2)
    S = -Z * (CosX * Sum2 - SinX * Sum1)
  ELSE
    C = 0: S = 0
  END IF
  IF X < 0 THEN C = -C: S = -S
END SUB

SUB GetKey (Choice, Maximum) STATIC
  A$ = ""
  WHILE A$ = "": A$ = INKEY$: WEND
  A = ASC(A$) - ASC("0")
  IF A > 0 AND A <= Maximum THEN
     Choice = A
  END IF
  PRINT Choice
END SUB

SUB GetNumber (Number) STATIC
   'if user entered RETURN without typing in a numeric value
   'then Number retains its original default value.
   INPUT Entry
   IF Entry <> 0 THEN Number = Entry
END SUB

SUB Intensity (X, Y, Z) STATIC
  Z = (X * X + Y * Y) / 2
END SUB

SUB LaserVector (X1, Y1, X2, Y2, ThePenWidth) STATIC
 'X1$,Y1$,X2$,Y2$ are local variables
 LaserXY X1, Y1, X1$, Y1$
 LaserXY X2, Y2, X2$, Y2$
 PRINT #Plotter, USING "PU; PW ##.##; PA"; ThePenWidth;
 PRINT #Plotter, X1$; Y1$; "; PD"; X2$; Y2$; "; PU;"
END SUB

SUB LaserXY (X, Y, X$, Y$) STATIC
   'Coordinate clipping
   IF X > DomainMax THEN
     X = DomainMax
   ELSEIF X < DomainMin THEN
     X = DomainMin
   END IF
   IF Y > Max THEN
     Y = Max
   ELSEIF Y < Min THEN
     Y = Min
   END IF
   X$ = STR$(INT(XBorderL + XScaleL * (X - DomainMin)))
   Y$ = STR$(INT(YBorderL + YScaleL * (Y - Min)))
END SUB

SUB Vector (X1, Y1, X2, Y2, TheColor) STATIC
 VideoXY X1, Y1, XX1, YY1
 VideoXY X2, Y2, XX2, YY2
 LINE (XX1, YY1)-(XX2, YY2), TheColor
END SUB

SUB VideoXY (X, Y, XX, YY) STATIC
   'Coordinate clipping
   IF X > DomainMax THEN
     X = DomainMax
   ELSEIF X < DomainMin THEN
     X = DomainMin
   END IF
   IF Y > Max THEN
     Y = Max
   ELSEIF Y < Min THEN
     Y = Min
   END IF
   XX = INT(XBorder + XScale * (X - DomainMin))
   YY = INT(YBorder + ScreenY - YScale * (Y - Min))
END SUB