/* Diffract.java

 * Computes the fresnel diffration pattern of various obstacles.

 * by John Talbot,  May 5,1997.

 */



class Table  // Optical table of components

{

  float a,b;

  Component c[3];



  void Table()

  {

    c[0]=new Beam();    // source component

    c[1]=new Wire();    // obstacle component

    c[2]=new Video();   // screen component

    float wavelength=(Source) c[0].wavelength;

    float source  =c[0].z;

    float obstacle=c[1].z;

    float screen  =c[2].z;



    float a1=obstacle-source;       //Source to Obstacle distance

    float a2=screen-obstacle;       //Obstacle to Screen distance

    float a3=screen-source;         //Source to Screen distance

    b=a3/a1;                        // ratio of big triangle to little triangle

    a=Math.sqrt(2/(wavelength*a2*b));   // Scaling down and SQR(2/(A3*L))

  }

  public static complex fresnel3(float origin, float x)

  {

    // a, b and parameters contain screen, obstacle and source dependence

    return(complex.fresnel(a*(b*origin-x)));

  }

  public shadowVector(float x1, float y1, float x2, float y2)

  {

     // project shadow of obstacle onto screen

     c[2].Vector(b*x1, b*y1, b*x2, b*y2);

  }

  public void add(component c)

  {

  }

  public void remove(component c)

  {

  }

}



class abstract Component

{

  float x=0, y=0, z=0;   // location of the component  (default = origin)

  String type;     // either 'Source' , 'Obstacle' or 'Screen'

  Component();

  {

    type="";  // required zero param constructor

  }

  Component(String type, float z)

  {

    this.z=z;   // default on the optical axis

    this.type=type;

  }

  setPosition(float x, float y, float z)

  {

    this.x=x;    // Can move the component

    this.y=y;

    this.z=z;

  }

  float getx()

  {

    return x;

  }

  float gety()

  {

    return y;

  }

  float getz()

  {

    return z;

  }

  String toString()

  {

    return type + " at ( " + x + ", " + y + ", " + z + ") ";

  }

}



class abstract Source extends Component

{

  float wavelength=550e-9;   // color  (default = 550 nm)

  String name="";            // either 'Spherical', 'Beam' or 'Plane'

  Source()

  {

    super("Source", 0);

  }

  Source(String name, float wavelength, float z)

  {

    super("Source", z);

    this.wavelength=wavelength;

    this.name=name;

  }

}



class abstract Obstacle extends Component

{

  String name="";    // either 'Edge', 'Slit', 'Wire', 'Two Slits' or 'Rectangle'

  Obstacle(String name, float z)

  {

    super("Obstacle", z);

    this.name=name;

  }

  // Set integration limits and compute intesity of diffraction patern

  abstract float amplitude(float x);  // obstacle dependent implementation

  // N.B.: should modify amplitude() in near future so it can return a

  // complex value for plotting the wavefunction phase



  abstract void shadow(Color EdgeColor);  // shadow cast by obstacle on screen



  String toString()

  {

    return name + super.toString();

  }

}



class abstract Screen extends Component

{

  float xMin=-20e-3;   // All Screens have these members in common

  float xMax= 20e-3;

  String name="";   // either 'Detector', 'Video' or 'VRML'

  Screen(String name, float z, float xMin, float xMax)

  {

    super("Screen", z);

    this.xMin=xMin;

    this.xMax=xMax;

    this.name=name;

  }

  abstract Vector(float x1, float y1, float x2, float y2);



}



// Linear detector scans across diffraction pattern and

// only produces a line graph

class Detector extends Screen

{

  int maxX=640;

  int maxY=480;

  float headRoom=2/100;   // 2 percent beyond extremums

  int xBorder=5;

  int yBorder=18;

  int screenX=maxX-2*xBorder;

  int screenY=maxY-2*yBorder;

  Detector()

  {

    super();

  }

  Dectector(float z, float xMin, float xMax)

  {

    super("Detector", z, xMin, xMax)

  }

}



// 2D Video detector scans a rectangular portion of the

// diffraction pattern and plots a color-coded intensity pixels

class Video extends Screen

{

  float yMin=-20e-3;   // extra dimension needed for 2D video detector

  float yMax= 20e-3;

  float a,b;

  Graphics g;          // Graphic context for output

  int maxX=320;

  int maxY=200;

  int maxPt=maxX;     // Increase this for super-sampling

  float aspect=1.3;

  float headRoom=0;    // 0 percent beyond extremums

  int xBorder=0;

  int yBorder=0;

  int screenX=maxX-2*xBorder;

  int screenY=maxY-2*yBorder;

  Color edgeColor=65;

  Color titleColor=64;

  Color borderColor=66;

  Color overlayColor=64;

  Color plotColor=63;

  // GreyScales();



  Video()

  {

    super("Video");

    g=getGraphics();

  }

  Video(float z, float xMin, float xMax, float yMin, float yMax)

  {

    super("Video", z, xMin, xMax);

    this.yMin=yMin;

    this.yMax=yMax;

    g=getGraphics();  // Need a Graphic context for plot

  }

  void Vector (float x1, float y1, float x2, float y2)

  {

    Integer xx1,yy1,xx2,yy2;

    // in case point is outside drawing area should handle

    // clipping and intersection with boundary here

    videoXY(x1,x1,xx1,yy1);

    videoXY(x2,x2,xx2,yy2);

    g.drawLine(java.lang.Integer.intValue(xx1),

               java.lang.Integer.intValue(yy1),

               java.lang.Integer.intValue(xx2),

               java.lang.Integer.intValue(yy2));

  }

  void videoXY (float x, float y, Integer xx, Integer yy)

  {  /* Map real world coordinates to screen coordinates */

    if(x>xMax) /* clipping */

      x=xMax;

    else if(x<xMin)

      x=xMin;

    if(y>yMax)

      y=yMax;

    else if(y<yMin)

      y=yMin;

    xx=xBorder+xScale*(x-xMin);

    yy=yBorder+screenY-yScale*(y-yMin);

  }

}



class Edge extends Obstacle

{

  float center=0;

  static complex c=new complex(0.5, 0.5);

  Edge()

  {

    super("Edge");

  }

  Edge(float z)

  {

    super("Edge", z);

  }

  complex amplitude(float x)

  {

    return complex.sub(c, Table.fresnel(center,x));

  }

  void shadow()

  {

    Table.shadowVector(xa,-java.lang.Float.MAX_VALUE,

                       xa,+java.lang.Float.MAX_VALUE);

  }

  String toString()

  {

    return super.toString() + " at position = " + center;

  }

}



class Slit extends Obstacle

{

  float width=1.3e-4;   // slit width

  float xa,xb;          // intermediate values to speed up computation

  Slit()

  {

    super("Slit");

    xb=width/2;

    xa=-xb;

  }

  Slit(float z)

  {

    super("Slit", z);

    xb=width/2;

    xa=-xb;

  }

  Slit(float z, float width)

  {

    super("Slit", z);

    this.width=width;

    xb=width/2;

    xa=-xb;

  }

  complex amplitude(float x)

  {

    complex z1=complex.fresnel(xa,x);

    complex z2=complex.fresnel(xb,x);

    return complex.Sub(z2,z1);

  }

  void shadow()

  {

    Table.shadowVector(xa,-java.lang.Float.MAX_VALUE,

                       xa,+java.lang.Float.MAX_VALUE);

    Table.shadowVector(xb,-java.lang.Float.MAX_VALUE,

                       xb,+java.lang.Float.MAX_VALUE);

  }

  String toString()

  {

    return super.toString() + " of width =" + width;

  }

}



class Wire extends Obstacle

{

  float width=1.3e-4;   // wire width

  float xa,xb;          // intermediate values to speed up computation

  Wire()

  {

    super("Wire");

    xb=width/2;

    xa=-xb;

  }

  Wire(float z)

  {

    super("Wire", z);

    xb=width/2;

    xa=-xb;

  }

  Wire(float z, float width)

  {

    super("Wire", z);

    this.width=width;

    xb=width/2;

    xa=-xb;

  }

  complex amplitude(float x)

  {

    Complez z=new complex(1.0, 1.0);

    complex z1=Table.fresnel(xa,x);

    complex z2=Table.fresnel(xb,x);

    return complex.sub(complex.add(z,z1),z2);

  }

  void shadow(Screen theScreen)

  {

    Table.shadowVector(xa,-java.lang.Float.MAX_VALUE,

                       xa,+java.lang.Float.MAX_VALUE);

    Table.shadowVector(xb,-java.lang.Float.MAX_VALUE,

                       xb,+java.lang.Float.MAX_VALUE);

  )

  String toString()

  {

    return super.toString() + " of width =" + width;

  }

}



class Rectangle extends Obstacle

{

  float width=1.3e-4;

  float height=2.5e-4;

  float xa,xb;          // intermediate values to speed up computation

  Rectangle()

  {

    super("Rectangle")

    xb=width/2;

    xa=-xb;

  }

  Rectangle(float z, float width, float height)

  {

    super("Rectangle", z);

    this.width=width;

    this.height=height;

    xb=width/2;

    xa=-xb;

  }

  complex amplitude(float x)

  {

    complex z1=Table.fresnel(xa,x);

    complex z2=Table.fresnel(xb,x);

    return complex.sub(z2,z1);

  }

  void shadow(Color edgeColor)

  {

    float x1,y1, x2,y2;

    x1=xa;

    x2=xb;

    y1=-height/2;

    y2= height/2;

    Table.shadowVector(x1,y1,x2,y1);

    Table.shadowVector(x2,y1,x2,y2);

    Table.shadowVector(x2,y2,x1,y2);

    Table.shadowVector(x1,y2,x1,y1);

  }

  String toString()

  {

    return super.toString() + " of width = " + width + " and height = " + height;

  }

}



class TwoSlits extends Obstacle

{ // Need GUI components to enter and display parameters

  float firstWidth=1e-4;    // width of slit 1

  float secondWidth=1e-4;   // width of slit 2

  float separation=2e-4;    // separation between the slits

  float xx,xz,xb,xa;        // intermediate values to speed up computation



  TwoSlits()

  {

    super("Two Slits");

    xx= separation/2;

    xz=-separation/2;

    xb= separation/2+secondWidth;

    xa=-separation/2-firstWidth;

  }

  TwoSlits(float z, float firstWidth, float secondWidth, float separation)

  {

    super("Two Slits", z);

    this.firstWidth=firstWidth;

    this.secondWidth=secondWidth;

    this.separation=separation;

    xx= separation/2;

    xz=-separation/2;

    xb= separation/2+secondWidth;

    xa=-separation/2-firstWidth;

  }

  complex amplitude(float x)

  {

      complex z1=Table.fresnel(xa,x);

      complex z2=Table.fresnel(xb,x);

      complex z3=Table.fresnel(xx,x);

      complex z4=Table.fresnel(xz,x);

      return complex.sub(complex.add(z2,z4),

                         complex.add(z1,z3));

  }

  void shadow()

  {

    Table.shadowVector(xa,-java.lang.Float.MAX_VALUE,

                       xa,+java.lang.Float.MAX_VALUE);

    Table.shadowVector(xx,-java.lang.Float.MAX_VALUE,

                       xx,+java.lang.Float.MAX_VALUE);

    Table.shadowVector(xz,-java.lang.Float.MAX_VALUE,

                       xz,+java.lang.Float.MAX_VALUE);

    Table.shadowVector(xb,-java.lang.Float.MAX_VALUE,

                       xb,+java.lang.Float.MAX_VALUE);

  }

  String toString()

  {

    return super.toString() + " with slit1 = " + firstWidth +

                              " wide, and slit2 = " + secondWidth +

                              " wide, separated by " + separation;

  }

}



class Spherical extends Source

{

  Spherical();

  {

    super("Spherical");

    wavelength=550e-9;

  }

  Spherical(float wavelenght, float z)

  {

    super("Spherical", z);

    this.wavelength=wavelength;

  }

}



class Beam extends Source

{

  float waist=1;    // angular beam spread  (default 1)

  Beam();

  {

    super("Beam")

    waist=1;

  }

  Beam(float wavelength, float z)

  {

    super("Beam", z);

  }

  Beam(float wavelength, float z, float waist)

  {

    super("Beam", z);

    this.waist=waist;

  }

  String toString()

  {

    return super.toString() + " and spread = " + waist;

  }

}



/*                    complex.java

 *

 * Common math routines for complex numbers and also computes specialized

 * complex functions such as the fresnel integrals.

 * by John Talbot,  May 5, 1997.

 */



class complex

{

  float real,imag;

  static float B1[12]=

  { 1.595769140,

   -0.000001702,

   -6.808568854,

   -0.000576361,

    6.920691902,

   -0.016898657,

   -3.050485660,

   -0.075752419,

    0.850663781,

   -0.025639041,

   -0.150230960,

    0.034404779};



  static float B2[12]=

  {-0.000000033,

    4.255387524,

   -0.000092810,

   -7.780020400,

   -0.009520895,

    5.075161298,

   -0.138341947,

   -1.363729124,

   -0.403349276,

    0.702222016,

   -0.216195929,

    0.019547031};



  static float B3[12]=

  { 0.000000000,

   -0.024933975,

    0.000003936,

    0.005770956,

    0.000689892,

   -0.009497136,

    0.011948809,

   -0.006748873,

    0.000246420,

    0.002102967,

   -0.001217930,

    0.000233939};



  static float B4[12]=

  { 0.199471140

    0.000000023

   -0.009351341

    0.000023006

    0.004851466

    0.001903218

   -0.017122914

    0.029064067

   -0.027928955

    0.016497308

   -0.005598515

    0.000838386};



  complex(float real, float imag)

  {

    this.real=real;

    this.imag=imag;

  }

  complex(float phase)

  {

     real=Math.cos(u);

     imag=Math.sin(u);

  }



  /* 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

   *

   * These integrals take a single real parameter and produce a single

   * complex number therefore this complex 'function' was placed in

   * the separate 'complex' class because it doesn't explicitely depend

   * on wavelength nor the distances to the screen, source and obstacle,

   * nor do they directly depend on the particular type of source,

   * obstacle or screen.

   */

  public static complex fresnel(float x)

  {

    complex c;

    float z;

    float u=Math.abs(x*x*Math.PI/2);

    float cosX=Math.cos(u);

    float sinX=Math.sin(u);

    if(u>4)

    {

      z=4/u;

      float sum3=0;

      float sum4=0;

      for(int i=11; i>0; i--)

      {

        sum3=(sum3+b3[i])*z;

        sum4=(sum4+b4[i])*z;

      }

      sum3=sum3+b3[0];

      sum4=sum4+b4[0];

      z=Math.sqrt(z);

      c=complex(0.5+z*(cosX*sum3+sinX*sum4),

                0.5-z*(cosX*sum4-sinX*sum3));

    }

    else if(u>0)

    {

      z=u/4;

      float sum1=0;

      float sum2=0;

      for(int i=11; i>0; i--)

      {

        sum1=(sum1+b1[i])*z;

        sum2=(sum2+b2[i])*z;

      }

      sum1+=b1[0];

      sum2+=b2[0];

      z=Math.sqrt(z);

      c=complex(z*(cosX*sum1+sinX*sum2),

               -z*(cosX*sum2-sinX*sum1));

    }

    else

    {

      c=complex(0,0);

    }

    if(x<0)

      return negate(c);

    else

      return c;

  }



  /* version of fresnel() using complex math functions

  public static complex fresnelb(float x)

  {

    complex c;

    float z;

    float u=Math.abs(x*x*Math.PI/2);

    complex phase=new complex(u);

    complex sum=new complex(0,0);

    if(u>4)

    {

      z=4/u;

      for(int i=11; i>0; i--)

        sum=scale(z, add(sum,c2[i]));

      sum.add(c2[i]);

      c=add(complex(0.5,0.5), scale(Math.sqrt(z),

                  multiply(phase,congugate(sum))));

    }

    else if(u>0)

    {

      z=u/4;

      for(int i=11; i>0; i--)

        sum=scale(z, add(sum,c1[i]));

      sum.add(c1[0]);

      c=scale(Math.sqrt(z), multiply(phase,congugate(sum)));

    }

    else

      c=complex(0,0);

    if(x<0)

      return negate(c);

    else

      return c;

  }

  */



  // *********** class methods *****************



  static complex add(complex x, complex y)

  {

    return complex(x.real+y.real,x.imag+y.imag);

  }

  static complex sub(complex x, complex y)

  {

    return complex(x.real-y.real,x.imag-y.imag);

  }

  static complex multiply(complex x, complex y)

  {

    return complex(x.real*y.real-x.

  }

  static complex intensity(complex x)

  {

    return(x.real*x.real+x.imag*x.imag);

  }

  static complex congugate(complex x)

  {

    return complex(x.real,-x.imag);

  }

  static complex negate(complex x)

  {

    return(-x.real,-x.imag);

  }

  static complex scale(float a, complex x)

  {

    return complex(a*x.real,a*x.imag);

  }



  //*********** instance methods ***************

  void set(float real, float imag)

  {

     this.real=real;

     this.imag=imag;

  }

  void add(complex x)

  {

    real+=x.real;

    imag+=x.imag;

  }

  void sub(complex x)

  {

    real-=x.real;

    imag-=x.imag;

  }

  void multiply(complex x, complex y)

  {

    real=real*x.real-imag*x.imag;

    imag=imag*x.real+real*x.imag;

  }

  float intensity()

  {

    return(real*real+imag*imag);

  }

  void congugate()

  {

    imag=-imag;

  }

  void negate()

  {

    real=-real;

    imag=-imag;

  }

  void scale(float a)

  {

    real=a*real;

    imag=a*imag;

  }



}



// ***************************************************************

// portions of old BASIC program partially converted



float Y[3300], PI, A, B, j[500], B1[11], B2[11], B3[11], B4[11];

float XBorderL, YBorderL, XScaleL, YScaleL, DomainMin, DomainMax, Min, Max;

float XBorder, YBorder, XScale, YScale, ScreenY, Plotter, MaxPt, XRange;



Constants();

// 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)



GaussWidth=-WG*WG/2;  //should compute this value

BesselZero=0.8*WG;    //rough initial estimate of first zero

OverlayFlag=1;        //if 1 then plot enveloppe overlay



Obstacle myObstacle=new Wire();  // can change the obstacle using the GUI



while(Plot!=Quit)

{

   MenuChoice();

   X=DomainMin;

   dX=XRange/(MaxPt-1);

   for(Pt=0; Pt<MaxPt; Pt++)  // compute and store results in an array

   {               // should be modified to plot it 'real-time'

     Y[Pt]=Convolve(myObstacle.amplitude(x),x);     //Function of X

     X+=dX;

   }



   if(Plot==Laser || Plot==Video)

     LaserPlot();

   elseif(Plot==TellAGraf)

     TellAGrafFile();

   elseif(Plot==Simulation)

   {

     if(Obstacle==Rectangle)

       Rectangular();

     else

       LaserSpot();

   }

}



'---------------------------------------------------------------------



public void MenuChoice()

{

 EdgeColor=10;

 TitleColor=14;

 BorderColor=5;

 OverlayColor=6;

 PlotColor=15



 Println("     Fresnel Diffraction          by John Talbot, Mar 16,1997");

 Println("     ------------------- ");

 Println("");

 Println(" Computer simulation of diffraction patterns produced by various ");

 Println(" simple obstacles using the Fresnel approximation.");

 Println(" (ref: Miles V. Klein 'Optics',1970 p.363-392, and AJP 40(1972)74)");

 Println("    1. None");

 Println("    2. WIRE");

 Println("    3. Edge");

 Println("    4. Slit");

 Println("    5. Two slits");

 Println("    6. Rectangular aperture");

 Println("    7. Quit");

 Println(" What type of obstacle ?");



 GetKey(Obstacle);

 if(Obstacle==Quit) exit(0);

 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);

 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 || Plot==TellAGraf)

  {

    Plotter++;

    Plotter$ = STR$(Plotter): FileName$ = CurrentDir$ + "Plot" + RIGHT$(Plotter$, LEN(Plotter$) - 1)

    System.out.println("  Will be saved to disk as % s", FileName);

    System.out.println();

  }

  if(Plot==Laser || Plot==Video) {

     PRINT "    1. Plot overlay of undisturbed laser intensity"

     PRINT "    2. No Overlay"

     PRINT " Plot Characteristics ?";

     GetKey OverlayFlag, 2#

     PRINT

  }



 else

  Plot==Simulation;

 }



 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=Math.sqrt(2*A1/(Wavelength*A2*(A1+A2)));  // Scaling down and sqrt(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)

 {

   System.out.println("Illegal Range ");

   exit(0);

 }



 switch(Plot)

 {

  case : Video

  case : Animation

   MaxPt=ScreenX;

   break;

  case : Laser

   MaxPt=ScreenXL;   //DO NOT not change this!

   break;

  case : TellAGraf

   MaxPt=ScreenXL;  //can be changed

   break;

 }



};



/*-------------------------------------------------------------*/

public void Constants()

{





  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 300DPI)



  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



  String XAxis = "POSITION (mm)";

  String YAxis = "INTENSITY (arbitrary units)";



  PlotPen=0.09;     // mm or  2 dots width

  TextPen=0.17;     // mm or  2 dots width

  TitlePen=0.34;    // mm or  4 dots width

  BorderPen=0.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;



  String CurrentDir="C:";  // Output directory

  Plotter=0;        // PCL5 and HPGL/2 plot file number for LaserJetIIIp.

}



void scale()

{

    max=-1E+20;

    min= 1E+20;

    for(pt=0; pt<maxPt; pt++)

    {

      y=y[Pt];

      if(y>max) max=y;

      if(y<min) min=y;

    }



    yRange=max-min;

    if(yRange==0) yRange=1;

    min=min-yRange*headRoom;

    max=max+yRange*headRoom;



    yRange=max-min;

    if(yRange==0) yRange=1;



    yScale=screenY/yRange;

    yScaleL=screenYL/yRange;



    xRange=domainMax-domainMin;

    xScale=screenX/XRange;

    xScaleL=screenXL/xRange;

}





/*-------------------------------------------------------------------*/

//Multiply the diffraction pattern by the Laser Beam Angular Envellope



public float Convolve(float intensity, float x)

{

   switch(Beam)

   {

     case: Gaussian

       return 10*intensity*Math.exp(x*x/GaussWidth);

     case: Bessel

       if(Plot!=Simulation) {

         Rho=3.832*x/BesselZero;

         BesselJ1(Rho, j1);

         j1=2*j1/Rho;

         return 10*intensity*j1*j1;

       }

       break;

     case: Uniform

       return 10*intensity;

    }

}



'-------------------------------------------------------------------------

'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



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





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



void Rectangular()

{

   Scale();

   MaxSpot=Max-YRange*HeadRoom   // reduce this for fainter detail



   Min=-Aspect*ScreenY/(2*XScale);

   Max=-Min;

   YRange=Max-Min;

   YScale=ScreenY/YRange;



   for(int j=0; j<ScreenY; j++)

   {

     Y=Aspect*(j-ScreenY/2)/XScale;     // in millimeters

     X=Y;

     Amplitude();

     YJ=256*YJ*Math.exp(Y*Y/GaussWidth)/MaxSpot;

     for(int i=0; i<ScreenX; i++)

     {

       X= INT(i*(MaxPt-1)/(ScreenX-1))

       Z=YJ*Y[X];

       if(Z>255) Z=255;  // set to maximum saturation

       PSET (i, j), Z;

     }

   }

}



'------------------------------------------------------

void LaserSpot()

{

   int X;

   Color Z;

   Graphics g = getGraphics();

   Scale();

   MaxSpot=Max-YRange*HeadRoom;  // reduce this for fainter detail

   if(Beam==Gaussian)

   {

     for(int j=0;j<ScreenY;j++)

     {

       Y=Aspect*(j-MAXY/2)/XScale;          // in millimeters

       YJ=64*Math.exp(Y*Y/GaussWidth)/MaxSpot;

       for(int i=0;i<ScreenX;i++)

       {

         X=i*(MaxPt-1)/(MAXX-1);

         Z=YJ*Y[X];

         if(Z>255) Z=255;     // set to maximum saturation

         g.setColor(Z);

         g.drawRect(i,j,1,1);

       }

     }

   }

   else if(Beam==Uniform)

     for(int i=0; i<ScreenX; i++)

     {

       X=i*(MaxPt-1)/(ScreenX-1);

       Z=256*Y[X]/MaxSpot;

       if(Z>255) Z=255;     // set to maximum saturation

       g.setColor(Z);

       g.drawLine(i,0,i,MAXY-1);

     }

   else if(Beam==Bessel)

   {

     for(int j=0; j<ScreenY;j++)

     {

       Y=Aspect*(j-ScreenY/2)/XScale;    // in millimeters

       for(int i=0; i<ScreenX; i++)

       {

         X=i*(MaxPt-1)/(ScreenX-1);

         XX=DomainMin+i*XRange/(ScreenX-1); // in millimeters

         Z=Y[X]/MaxSpot;

         Rho=3.832*Math.sqrt(XX*XX+Y*Y)/BesselZero;

         BesselJ1(Rho,j1);

         j1=2*j1/Rho;

         Z=256*Z*j1*j1;

         if(Z>255) Z=255;     // set to maximum saturation

         g.setColor(Z);

         g.drawRect(i,j,1,1);

       }

     }

   }



   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



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 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



public float Intensity(float X, float Y)

{

  return (X*X+Y*Y)/2;

};



/* Should output this to a VRML file format in future */

LaserVector(float X1, float Y1, float X2, float Y2, float ThePenWidth)

{

 String X1s,Y1s,X2s,Y2s;

 LaserXY(X1,Y1,X1s,Y1s);

 LaserXY(X2,Y2,X2s,Y2s);

 System.out.println("PU; PW %3.2f ; PA %s %s ; PD %s %s ; PU;",

                     ThePenWidth,X1s,Y1s,X2s,Y2s);

}



void LaserXY (float X, float Y, String X, String Y)

{  /* Map real world coordinates to screen coordinates */

   if(X>DomainMax) /* clipping */

     X=DomainMax;

   else if(X<DomainMin)

     X=DomainMin;

   if(Y>Max)

     Y=Max;

   else if(Y<Min)

     Y=Min;

   String X=new( (int) XBorderL+XScaleL*(X-DomainMin) )

   String Y=new( (int) YBorderL+YScaleL*(Y-Min) )

}