'---------------------------------------------------------------
'                   ANIMATE.BAS
'
'  http://www.achilles.net/~jtalbot/amateur/diffract
'
'  by John Talbot, jtalbot@achilles.net  Apr 19,1997.
'
'  This version works only on an amiga computer.
'
' Computes the fresnel diffration pattern of various obstacles
' and plots them to the screen as graphs and/or intensity coded
' colors. There is an option to (parts of this program are from
' SUPERGRAPH a rudimentary function plotter with many options).
'
'----------------------------------------------------------------

Main:
DEFDBL A-Z  
GOSUB SetConstants

CHDIR "ram:bmaps"
LIBRARY "graphics.library"    'PolyDraw&
CHDIR "df0:"
RANDOMIZE TIMER

FlagOn Mode,FASTPLAYBACK
MaxPt=150 : MaxPt2=18 : GOSUB InitialiseSuperGraph   
GOSUB Compute

END
'------------------------------------------
'  Program specific information
'---------------------------------------------------       
Compute:

  DIM SHARED A,B,B1(11),B2(11),B3(11),B4(11)
  GOSUB Coefficients

  None   = 0       'Type of obstacle
  Wire   = 1      
  Edge   = 2
  Slit   = 3
  TwoSlit= 4 

  Uniform = 0      'Types of Beam
  Gaussian = 1
  Bessel = 2
                   'All distances are in millimeter (mm)
  X0 = 20          'Scan sreen position from -X0 to +X0 mm
                   'AJP 40(1972)74           -30 to +30 mm
                   'Actual guassian outline  -80 to +80 mm
  Obstacle=Wire    'Type of aperture 
  Beam=Gaussian    'Type of laser beam
  '----------------------------------------
  Mode=NORMAL OR THIN     'LOGX 'OR LOGY ' or make Hires=1
  Y$="Amplitude"

  'First Parameter (plotted as X-Axis)
  X$="Position"
  DomainMin=-X0
  DomainMax=+X0   : IF DomainMin>DomainMax THEN END

  'Second Parameter (animation)
  DomainMin2=.1
  DomainMax2=.9   : IF DomainMin2>DomainMax2 THEN END

Restart:
  GOSUB ResetX2      'reset X2,dx2,LogDx2
  LOCATE 1,1:PRINT Pt2
  FOR Pt2=0 TO MaxPt2-1
    LOCATE 1,1:PRINT MaxPt2-Pt2
    GOSUB TheConstants   'using X2 we could Sweep A2, A1, XD
    GOSUB ResetX         'reset X,dx,LogDx
    FOR Pt=0 TO MaxPt-1    
      XJ=X : GOSUB Amplitude
      Y(Pt2,   Pt) = YJ    'Function of X
      Y(MaxPt2,Pt) = X     'XCoordinate    
      X=X+dx : dx=dx*LogDx    
    NEXT Pt
    Y(Pt2,MaxPt)=X2     'save X2
    ScaleGroup(Pt2)=0   'Pt2 scale each graph individually
    Y(Pt2,MaxPt)=X2  
    X2=X2+dX2 : dX2=dX2*LogDx2   
  NEXT Pt2
 '----------------------

  
  'Now using SuperGRAPH routines  
  FlagOn Mode,FASTPLAYBACK   'like an oscilloscope !!
  GOSUB GetScale
  CLS
  FlagOn Mode,GRID
  VTic=0
  'VTics 0#,VTic,PI   
  HTic=0
  HTics 0#,HTic,Max(0)
  HAxis 0#,Tol,3#
  HAxis 0#,HighTol,3#
  GOSUB PlotGraph

  END
  FlagOff Mode,FASTPLAYBACK
  GOTO Restart
RETURN    

'---------------------------------------------------------------------
TheConstants:
                   'all units are in mm (millimeters)

  WG = 52.1        'Gaussian envellope width at intensity = 13.1%
  A1 =   88 '+X2   'Source to Obstacle distance
  A2 = 1469 '-X2   'Obstacle to Screen distance
  A3 = 6.3282E-04  'Laser wavelength;  632.82 nm (monomode only)
  XD = X2 '.13     'If Obstacle = Wire    then XD = Diameter  
                   '            = Edge    then XD = Center
                   '            = Slit    then XD = Width
                   '            = TwoSlit then XD = Separation
  XA = 0           'IF Obstacle = TwoSlit then XA = First  slit width
  XC = 0           'IF Obstacle = TwoSlit then XC = Second slit width

  A  = SQR(2*A1/(A3*A2*(A1+A2)))      'Scaling down and SQR(2/(A3*L))
  B  = (A1+A2)/A1                     'Scaling up 
  IF Obstacle=TwoSlit THEN
    XX = XD/2    :  XZ =-XX
    XB = XX+XC   :  XA =-(XX+XA)
  ELSE
    XB = XD/2
    XA =-XB
  END IF
  GaussWidth = -WG * WG / 2
  BesselZero = WG / 4       'rough initial estimate of first zero

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 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
      Intensity 1, 0, YJ
    END IF

    IF Beam = Gaussian THEN
       YJ = YJ * EXP(X * X / GaussWidth)
    ELSEIF Beam = Uniform THEN
       YJ = 1 * YJ
    END IF
   
  RETURN
END
'-----------------------------------------------------------------
'Use the correct screen to aperture projection factors

SUB Fresnel(C,S,XOrigin,XPos) STATIC
  X = A*(B*XOrigin-XPos)
  Fresnel3 C,S,X
END SUB
'-----------------------------------------------------------------
'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 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
'---------------------------------------------------------
Coefficients:

  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 Coefficients
  FOR i=0 TO 11
    READ B1(i),B2(i),B3(i),B4(i)
  NEXT i
RETURN

'--------------------------------------------------------------
'Compute then fresnel diffration intensity of the complex amplitude

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

'-------END of COMPUTE part-----------------------------------


'----------------SuperGRAPH----------------------------
'------------------------------------------------------
'(To be merged to a specific application)
'  written by John Talbot ;  Dec87  ;Apr88  ;Jul 2,1990
'  Reduced in Complexity Feb 7,1993
'--------------------------------------------
ResetX:
  IF DomainMin>DomainMax THEN t=DomainMax : DomainMax=DomainMin:DomainMin=t
  XRange =DomainMax -DomainMin
  IF XRange =0 OR MaxPt =1 THEN PRINT "Domain Error":END
  IF (DomainMax<=0 OR DomainMin<=0) AND (Mode AND LOGX)<>0 THEN 
     FlagOff Mode,LOGX
     PRINT "LogX Mode Canceled :";DomainMin;DomainMax
  END IF
  IF Mode AND LOGX THEN 
    dLx=(LOG(DomainMax)-LOG(DomainMin))/(MaxPt-1) 
    dx =DomainMin*(EXP(dLx)-1)
  ELSE 
    dLx=0
    dx= XRange/(MaxPt-1)
  END IF
  LogDx=EXP(dLx)
  X=DomainMin

RETURN

ResetX2:
  IF DomainMin2>DomainMax2 THEN t=DomainMax2 : DomainMax2=DomainMin2:DomainMin2=t
  XRange2=DomainMax2-DomainMin2
  Intervals=MaxPt2-1
  IF Intervals>0 THEN
     IF XRange2=0 THEN PRINT "Warning X2 range is zero":END 
     IF (DomainMax2<=0 OR DomainMin2<=0) AND (Mode AND LOGX2)<>0 THEN 
        FlagOff Mode,LOGX2
        PRINT "LogX2 Mode Canceled :";DomainMin2;DomainMax2
     END IF
     IF Mode AND LOGX2 THEN 
       dLx2=(LOG(DomainMax2)-LOG(DomainMin2))/Intervals 
       dX2 =DomainMin2*(EXP(dLx2)-1)
     ELSE    'if only one AND in LOGX2 Mode
       dLx2=0
       dX2= XRange2/Intervals
     END IF
  ELSE       'if only one X2
     dLx2=0
     dX2=0
  END IF
  LogDx2=EXP(dLx2)
  X2=DomainMin2
  'To Prevent Plotting until you have called Getscale
RETURN

'---------------------------------------------
SetConstants:

DIM SHARED MaxPt2,PI
LARGE=1E+20     'Very Large number
PI=3.141592653589# 

'SUPERGRAPH Modes
DIM SHARED Status,Mode,LOGX,LOGY,THIN,GRID
NORMAL = 0*2^0
LOGX   = 1*2^3   ' ex: LOGX , LOGX|LOGY , LOGY 
LOGY   = 1*2^4
LOGX2  = 1*2^5
GRID   = 1*2^6

THIN   = 1*2^7   'Vector thickness
DOT    = 1*2^8   'Points are emphasized

HISTORY= 1*2^14     'this is for FastPlot Option
FASTPLAYBACK=1*2^13 'this turns on FastPlot

'SUPERGRAPH Status signals
FASTGRAPH =1*2^1
FASTALLOC =1*2^2
FIRSTCALL =1*2^3

FlagOff Status,FASTALLOC
FlagOn  Status,FIRSTCALL

RETURN
'---------------------------------------------

ResetSuperGraph:   'If NEW MaxPt,MaxPt2
  ERASE Y,Min,Max,YRange,YScale,X,ScaleGroup,MyColor
  IF Status AND FASTALLOC THEN
    ERASE YFast%
    FlagOff Status,FASTALLOC   'only place to turn it off
  END IF
RETURN

InitialiseSuperGraph:

IF (Status AND FIRSTCALL)=0 THEN GOSUB ResetSuperGraph
IF MaxPt<3 THEN MaxPt=3
IF MaxPt2<1 THEN MaxPt2=1
DIM SHARED X$,Y$,DomainMin,DomainMax,XScale,XRange
DIM SHARED ScreenX,ScreenY,Xborder,YBorder,Aspect
IF Mode AND FASTPLAYBACK THEN FlagOn Status,FASTALLOC

IF Status AND FASTALLOC THEN 
  MaxFrame=MaxPt2-1    
  DIM YFast%(2*MaxPt+2,MaxFrame)
END IF
DIM SHARED Y(MaxPt2,MaxPt)  'maximum screen points
DIM SHARED Min(MaxPt2),Max(MaxPt2)
DIM SHARED YRange(MaxPt2),YScale(MaxPt2)

DIM SHARED X(MaxPt)       'not necessary
DIM SHARED ScaleGroup(MaxPt2)

MAXX=610       'Maximum Available Plotting Area
IF Hires=1 THEN
 MAXY=385
 SCREEN 1,640,400,2,4
 WINDOW 1,"High Resolution",(0,0)-(MAXX,MAXY),15,1
 WINDOW OUTPUT 1
ELSE
 MAXY=190
END IF

Aspect=.44
PLANES=2     
'-----------Set Up Colors-----------------------
MaxColors=2^PLANES    'This Should Be a Parameter
DIM SHARED BLACK,YELLOW,BLUE,GREEN
DIM MyColor(MaxColors)
IF MaxColors=4 THEN
 BLACK=0:YELLOW=1 :BLUE=2 :GREEN=3 
 PALETTE BLACK  ,0 , 0, 0
 PALETTE YELLOW ,1 , 1, 0
 PALETTE BLUE   ,0 ,.3, 1
 PALETTE GREEN  ,0 , 1, .3

 MyColor(0)=GREEN
 MyColor(1)=BLUE
 MyColor(2)=GREEN
 MyColor(3)=YELLOW
ELSE
'  Note: if called again it will close and reopen screen 
' SCREEN (MAXX,MAXY)
' WINDOW
 BLACK=0:YELLOW=1 :BLUE=2 :GREEN=3 'defined in color sub
 FOR C=0 TO MaxColors-1
   MyColor(C)=C
 NEXT C
 MyColor(0)=GREEN
END IF
RastPort&=WINDOW(8)

FlagOff Status,FASTGRAPH
FlagOff Status,FIRSTCALL

RETURN
'-------FLAG MANIPULATORS----------------------

SUB FlagOn (StatusReg,Bit) STATIC
  StatusReg= StatusReg OR Bit  
END SUB

SUB FlagOff (StatusReg,Bit) STATIC
  StatusReg=(StatusReg OR Bit) XOR Bit  
END SUB

'-----------SUBROUTINES-------------------------
PlotGraph:
'Draw Graph Perimeter
  x1=DomainMin:y1=Min(0) : MapXY  0#,x1,y1
  X2=DomainMax:y2=Max(0) : MapXY  0#,X2,y2
  LINE (x1-1,y1+1)-(X2+1,y2-1),1,B   'make graph Box 
  LINE (x1-2,y1+1)-(X2+2,y2-1),1,B  

IF Status AND FASTGRAPH THEN  'we have Valid YFast%
  FOR Frame%=0 TO MaxFrame
    col=MyColor(Frame% MOD MaxColors)
    COLOR col,BLACK
    GOSUB VectorPlot    
  NEXT Frame%
ELSE    'if no we must compute every point  
  FOR Pt2=0 TO MaxPt2-1
     GOSUB DoGraph
  NEXT Pt2
  IF Mode AND FASTPLAYBACK THEN FlagOn Status,FASTGRAPH  'we now have YFast%      
END IF

IF Mode AND FASTPLAYBACK THEN GOSUB PlayBackFast

RETURN

DoGraph:
     Frame=Pt2

     LastY=Y(Pt2,0) :LastX=Y(MaxPt2,0) 
     MapXY Pt2,LastX,LastY
    
     col=MyColor(Frame MOD MaxColors)

     IF Mode AND FASTPLAYBACK THEN
       YFast%(0,Frame)=LastX 
       YFast%(1,Frame)=LastY
     END IF
     
     IF Mode AND DOT THEN   'Dot on Experimental Points 
       LINE (LastX-1,LastY-1)-STEP(2,2),col,bf    
     END IF

    FOR Pt=1 TO MaxPt-1   'start at next point (one)

      Y=Y(Pt2,Pt)  :X=Y(MaxPt2,Pt)  
      MapXY Pt2,X,Y    

      LINE (LastX,LastY)-(X,Y),col    
      IF Mode AND FASTPLAYBACK THEN
          YFast%(2*Pt  ,Frame)=X 
          YFast%(2*Pt+1,Frame)=Y
      END IF
      IF (Mode AND THIN)=0 THEN    'double line (in 640)
         LINE (LastX+1,LastY)-(X+1,Y),col
      END IF
      IF Mode AND DOT THEN   'Dot on Experimental Points 
        LINE (X-1,Y-1)-STEP(2,2),col,bf    
      END IF
      LastX=X:LastY=Y
    NEXT Pt

RETURN

'------Get Scale for Graphs and SetMode ( NORMAL ...)
GetScale:
  IF ((Mode AND FASTPLAYBACK)<>0 AND (Status AND FASTALLOC)=0) THEN
     PRINT "Sorry! FASTPLAYBACK Not Available,You Must set "
     PRINT "Mode=FASTPLAYBACK before calling InitialiseSuperGraph"
     PRINT "(This is Optional since YFast%(,) take lots of memory!)"
     FlagOff Mode,FASTPLAYBACK
  END IF
  'reset FASTGRAPH since we must now recompute YFast%()
  FlagOff Status,FASTGRAPH
  GOSUB SetMode
  GOSUB Scale
RETURN

'---------Scale Graph according to Extremum (X and Y)------
Scale:
  HeadRoom=3/100   ' 3% beyond extremums

  FOR Pt2=0 TO MaxPt2-1
    Max=-LARGE
    Min=+LARGE        
    FOR Pt=0 TO MaxPt-1
      Y=Y(Pt2,Pt)
      IF Y>Max THEN Max=Y
      IF Y<Min THEN Min=Y
    NEXT Pt
    YRange=Max-Min 
      
    IF YRange=0 THEN           'making shure YRange<>0
        IF (Max=0 AND Min=0) THEN
          YRange=1
        ELSEIF Max=Min THEN
          YRange=.99*Max     '99% so Min>0
        END IF
        Max= Max + YRange
        Min= Min - YRange
    ELSE    '10% Buffer on Extremums
        Max= Max + YRange*HeadRoom 
        IF Mode AND LOGY THEN
         'if LOGY and min>0 : HeadRoom ONLY if after min>0 
           IF NOT(Min>0 AND (Min-YRange*HeadRoom)<=0) THEN
             Min= Min - YRange*HeadRoom  
           END IF
        ELSE
             Min= Min - YRange*HeadRoom  
        END IF
    END IF    

    IF (Min<=0 OR Max<=0) AND (Mode AND LOGY)<>0 THEN
      FlagOff Mode,LOGY   
      PRINT "Canceling LOGY Mode : Negative or zero values."
      PRINT "At Point:";Pt2;Min;Max
    END IF
    
    Max(Pt2)=Max
    Min(Pt2)=Min
  NEXT Pt2


  '/* IMPORTANT: If DataSets are on Same Scale You MUST */
  '/* set ALL Min() and Max() to common Extremums       */

  'Select commom Y Scales if necessary 
  FOR i=0 TO MaxPt2-1
   FOR j=0 TO MaxPt2-1
      IF ScaleGroup(i)=ScaleGroup(j) THEN
         IF Max(i)>Max(j) THEN Max(j)=Max(i)
         IF Min(i)<Min(j) THEN Min(j)=Min(i)
      END IF
   NEXT j 
  NEXT i

  'Set Correct Plotting Scales      
  FOR Pt2=0 TO MaxPt2-1   
    YRange(Pt2)=Max(Pt2)-Min(Pt2)
    IF Mode AND LOGY THEN  
      Y=LOG(Max(Pt2))-LOG(Min(Pt2))
    ELSE
      Y=YRange(Pt2)
    END IF
    YScale(Pt2)=ScreenY/Y
  NEXT Pt2
  '.................................
  'Now take care of X-Parameters

    Max=-LARGE
    Min=+LARGE        
    FOR Pt=0 TO MaxPt-1
      Y=Y(MaxPt2,Pt)
      IF Y>Max THEN Max=Y
      IF Y<Min THEN Min=Y
    NEXT Pt
  Max(MaxPt2)=Max
  Min(MaxPt2)=Min

  XRange=Max(MaxPt2)-Min(MaxPt2)
  DomainMin=Min(MaxPt2)
  DomainMax=Max(MaxPt2)
    
  'Mode=LOGX is illegal if any X Values are 0 or negative
  IF DomainMin<=0 THEN Mode=(Mode OR LOGX) XOR LOGX

  IF Mode AND LOGX THEN
   XScale=ScreenX/( LOG(DomainMax)-LOG(DomainMin) )
  ELSE
   XScale=ScreenX/XRange     'same scale for all DataSets
  END IF
RETURN
'-------------MODE INITIALISATION-------------------
SetMode:
   Xborder=5  '50  '80   'should be a parameter
   YBorder=9  '10  '20
   IF 2*Xborder>MAXX OR 2*YBorder>MAXY THEN END
   ScreenX=MAXX-2*Xborder    'effective Plotting area
   ScreenY=MAXY-2*YBorder
RETURN

'---------PLOT Graph Major Axes-------------------
SUB VAxis (DataSet,HPos,MyColor) STATIC
  'LOCALS: X1,Y1,X2,Y2
 IF DataSet<MaxPt2 THEN
    IF DomainMin<HPos AND HPos<DomainMax THEN  
      x1=HPos : y1=Min(DataSet): MapXY DataSet,x1,y1
      X2=HPos : y2=Max(DataSet): MapXY DataSet,X2,y2
      LINE (x1,y1)-(X2+1,y2),MyColor,bf
    END IF
 END IF
END SUB

SUB HAxis (DataSet,VPos,MyColor) STATIC
  'LOCALS: X1,Y1,X2,Y2
 IF DataSet<MaxPt2 AND Min(DataSet)<VPos AND VPos<Max(DataSet) THEN
     x1=DomainMin : y1=VPos : MapXY DataSet,x1,y1
     X2=DomainMax : y2=VPos : MapXY DataSet,X2,y2
     LINE (x1,y1)-(X2,y2),MyColor
 END IF
END SUB
                                                     
'------AXIS 'TICKING' or gradations---------------------
'INPUT: DataSet,TicSize (0 for AUTO),XPos/YPos,Grid(BOOLEAN)

SUB VTics (DataSet,TicSize,XPos) STATIC
  'LOCALS: k,X,Y

 IF DataSet<MaxPt2 THEN
   IF TicSize=0 THEN FindTicSize YRange(DataSet),TicSize
   k=INT(Min(DataSet)/TicSize)+1    ' also good for Min<0
   WHILE k*TicSize<Max(DataSet)
      IF Mode AND GRID THEN 
         HAxis DataSet,k*TicSize,BLUE
      ELSE
         X=XPos :Y=k*TicSize :MapXY DataSet,X,Y        
         LINE (X,Y)-STEP(6,0),GREEN,bf           
      END IF
      k=k+1
   WEND          
 END IF

END SUB              

SUB HTics (DataSet,TicSize,YPos) STATIC
  'LOCALS: k,X,Y,Min,Max,Range

 IF DataSet<MaxPt2 THEN
   Min=DomainMin
   Max=DomainMax
   Range=Max-Min
   IF TicSize=0 THEN FindTicSize Range,TicSize
   k=INT(Min/TicSize)+1    
   WHILE k*TicSize<Max
      IF Mode AND GRID THEN     
         VAxis DataSet,k*TicSize,BLUE      
      ELSE   'We must generate Tics of 4 pixels long
         mag=+1/4
         x1=k*TicSize : y1=YPos   : MapXY DataSet,x1,y1        
         X2=k*TicSize : y2=Min(0)      
         IF YPos=Min(0) THEN y2=Max(0): mag=-mag
         MapXY DataSet,X2,y2        
         X2=(X2-x1): y2=(y2-y1)
         mag=mag*SQR(X2*X2+y2*y2)
         X2=X2/mag : y2=y2/mag  
         LINE (x1  ,y1)-STEP(X2  ,y2),YELLOW
         IF (Mode AND THIN)=0 THEN
            LINE (x1+1,y1)-STEP(X2,y2),YELLOW
         END IF
      END IF
      k=k+1
   WEND          

 END IF

END SUB

SUB FindTicSize (Range,TicSize) STATIC
                'INPUT,OUTPUT 
    'LOCALS: Temp,Mantissa,Exponent
    TicSize=.5    'default  
    temp=LOG(Range)/LOG(10)      'note: Range is > 0
    Exponent=INT(temp)
    Mantissa=10^(temp-Exponent)
    IF Mantissa>=1!  AND Mantissa<2.5  THEN TicSize=.2   
    IF Mantissa>=2.5 AND Mantissa<5!   THEN TicSize=.5 
    IF Mantissa>=5!  AND Mantissa<9.99 THEN TicSize=1 
    TicSize=TicSize*10^Exponent
END SUB

'-------Scaled Screen OUTPUT Mapping-------------------
A:

SUB MapXY (DataSet,X,Y) STATIC
 IF Mode AND LOGX THEN 
  X=Xborder + XScale*(LOG(X)-LOG(DomainMin))
 ELSE
  X=Xborder + XScale*(X-DomainMin)
 END IF
 IF Mode AND LOGY THEN
   Y=YBorder+ScreenY-YScale(DataSet)*(LOG(Y)-LOG(Min(DataSet)) )
 ELSE
   Y=YBorder+ScreenY-YScale(DataSet)*(Y-Min(DataSet))   
 END IF
END SUB

'------------fast playback (animate a sequence of graphs)----------
PlayBackFast:

COLOR YELLOW,BLACK
FlagOff Mode,HISTORY
MaxDelay=50 : FrameStep=1 : MaxFrame=MaxPt2-1 : Frame%=0
A$="" : Reverse=+1    'if Neg then bounceback else restart
Monitor$="Source 0 mm  Obstacle ####.# mm  Dectector ####.# mm  Width #.###-#.### mm"
PRINT USING,Monitor$;A1;A1+A2;DomainMin2;DomainMax2;
WHILE INKEY$="" : WEND

Monitor$="Source 0 mm  Obstacle ####.# mm  Dectector ####.# mm  Width #.### mm"
WHILE( A$<>"q" )
 A$=""
 WHILE (A$="") 
     A$=INKEY$ : GOSUB NoClean : GOSUB VectorPlot
     LOCATE 1,1: PRINT USING Monitor$;A1;A1+A2;DomainMin2+dX2*Frame%;
     FOR Delay=0 TO MaxDelay : NEXT Delay
     Frame%=Frame%+FrameStep
     IF Frame%<0 THEN 
       'There is a bug here which prevents a full sweep of Frame%
       IF Reverse<0 THEN Frame%=1 ELSE Frame%=MaxFrame
       FrameStep=Reverse*FrameStep
       GOSUB Clean
     ELSEIF Frame%>MaxFrame THEN 
       IF Reverse<0 THEN Frame%=MaxFrame ELSE Frame%=0
       FrameStep=Reverse*FrameStep
       GOSUB Clean
     END IF
 WEND
 IF A$="r" THEN 
    Reverse=-1*Reverse
 ELSEIF A$=" " THEN 
    WHILE INKEY$="" : WEND
 ELSEIF A$="t" THEN 
    LOCATE 1,2: INPUT "What Rate";MaxDelay
    MaxDelay=ABS(MaxDelay)+1
 ELSEIF A$="h" THEN
    Mode=Mode XOR HISTORY
 END IF
WEND 
RETURN
'----------VECTOR-PLOTTER-------------------
VectorPlot:
  CALL Move&(RastPort&,YFast%(0,Frame%),YFast%(1,Frame%))
  CALL PolyDraw&(RastPort&,CINT(MaxPt-1),VARPTR(YFast%(2,Frame%))   )  
  'PSET (YFast%(0,Frame%),YFast%(1,Frame%)),1
  'Index%=2
  'FOR Pt=1 TO MaxPt-1   'start at next point (one)
  '   LINE -(YFast%(Index%,Frame%)-(Index%+1,Frame%),1
  '   Index%=Index%+2
  'NEXT Pt
RETURN

NoClean:
  IF (Mode AND HISTORY)=0 THEN
    LINE (Xborder,YBorder+1)-STEP(ScreenX-1,ScreenY-2),0,bf 
  END IF
RETURN   
'--------------Clear window----------
Clean:
  IF Mode AND HISTORY THEN
    LINE (Xborder,YBorder+1)-STEP(ScreenX-1,ScreenY-2),0,bf
  END IF
RETURN