DECLARE EXTERNAL

EXTERNAL DivisionIndex%:(XDivisionNum%,YDivisionNum%,NumXDivisions%,NumYDivisions%)
EXTERNAL DrawLine:
EXTERNAL F:(X,Y)

PROC Plot3D:
  GLOBAL BoxX1%
  GLOBAL BoxX2%
  GLOBAL BoxY1%
  GLOBAL InterceptCountMod2%
  LOCAL AspectRatio
  LOCAL Bias
  LOCAL BoxDeltaX
  LOCAL BoxDeltaY
  LOCAL BoxNum1%
  LOCAL BoxNum2%
  LOCAL BoxX%(4)
  LOCAL BoxXIntercept
  LOCAL BoxY%(4)
  LOCAL BoxYMax%
  LOCAL BoxYMin%
  LOCAL BoxYOffset
  LOCAL Color%(1000)
  LOCAL ColorMax%
  LOCAL ColorMin%
  LOCAL ColorNum%
  LOCAL CosRotation
  LOCAL CosTilt
  LOCAL DarkestGray%
  LOCAL DeltaX
  LOCAL DeltaY
  LOCAL DeltaZ
  LOCAL DivIndex%
  LOCAL Finished%
  LOCAL Fraction
  LOCAL idPlot%
  LOCAL KeyIndex1%
  LOCAL KeyIndex2%
  LOCAL Largest
  LOCAL LightestGray%
  LOCAL LightX
  LOCAL LightY
  LOCAL LightZ
  LOCAL Magnitude
  LOCAL MaxYOut%
  LOCAL MaxZOut%
  LOCAL NormalX
  LOCAL NormalY
  LOCAL NormalZ
  LOCAL NumXDivisions%
  LOCAL NumXPrimes%
  LOCAL NumYDivisions%
  LOCAL PixelsPerUnit
  LOCAL Response$(1)
  LOCAL Radians
  LOCAL RadiansPerDegree
  LOCAL Rotation
  LOCAL SinRotation
  LOCAL SinTilt
  LOCAL Smallest
  LOCAL SortLeft%
  LOCAL SortRight%
  LOCAL SortT1
  LOCAL SortT2%
  LOCAL SortT3%
  LOCAL Tilt
  LOCAL VertexX(4)
  LOCAL VertexY(4)
  LOCAL VertexZ(4)
  LOCAL X
  LOCAL XDivisionIndex%(1000)
  LOCAL XDivisionNum%
  LOCAL XEye
  LOCAL XMax
  LOCAL XMin
  LOCAL XPrime(1000)
  LOCAL XPrimeMax
  LOCAL XPrimeNum%
  LOCAL XRotated
  LOCAL Y
  LOCAL YCenter
  LOCAL YDivisionIndex%(1000)
  LOCAL YDivisionNum%
  LOCAL YMax
  LOCAL YMin
  LOCAL YOffset
  LOCAL YOutMax
  LOCAL YPrime(1000)
  LOCAL YPrimeMax
  LOCAL YPrimeMin
  LOCAL Z
  LOCAL ZCenter
  LOCAL ZOffset
  LOCAL ZOutMax
  LOCAL ZPrime(1000)
  LOCAL ZPrimeMax
  LOCAL ZPrimeMin
  CLS
  MaxYOut%=639
  MaxZOut%=239
  PRINT "                             Three Dimensional Plot"
  PRINT
  PRINT
  PRINT
  PRINT "Smallest value for x? ";
  INPUT XMin
  PRINT "Largest value for x? ";
  INPUT XMax
  PRINT "Smallest value for y? ";
  INPUT YMin
  PRINT "Largest value for y? ";
  INPUT YMax
  DO
    PRINT "Number of divisions for x? ";
    INPUT NumXDivisions%
    IF NumXDivisions% <= 1 OR NumXDivisions% > 1000
      PRINT "? there must be at least 2 and no more than 1000 divisions"
    ENDIF
  UNTIL NumXDivisions% > 1 AND NumXDivisions% <= 1000
  DO
    PRINT "Number of divisions for y? ";
    INPUT NumYDivisions%
    IF NumYDivisions% <= 1 OR NumYDivisions% > 1000/NumXDivisions%
      PRINT "? there must be at least 2 and no more than ";1000/NumXDivisions%;" divisions"
    ENDIF
  UNTIL NumYDivisions% > 1 AND NumYDivisions% <= 1000/NumXDivisions%
  DO
    PRINT "Rotation about the z-axis (degrees)? ";
    INPUT Rotation
    PRINT "Tilt about the resulting y-axis (degrees)? ";
    INPUT Tilt
    PRINT "After the plot is displayed, press a key to continue."
    PRINT "Evaluating function..."
    RadiansPerDegree=ATAN(1.0)/45.0
    Radians=Tilt*RadiansPerDegree
    CosTilt=COS(Radians)
    SinTilt=SIN(Radians)
    Rotation=-(Rotation+90)
    Radians=Rotation*RadiansPerDegree
    CosRotation=COS(Radians)
    SinRotation=SIN(Radians)
    X=XMin
    Y=YMin
    Z=F:(X,Y)
    XRotated=XMin*CosRotation+YMin*SinRotation
    YPrimeMin=-XMin*SinRotation+YMin*CosRotation
    ZPrimeMin=-XRotated*SinTilt+Z*CosTilt
    YPrimeMax=YPrimeMin
    ZPrimeMax=ZPrimeMin
    XPrimeMax=XRotated*CosTilt+Z*SinTilt
    DeltaX=NumXDivisions%
    DeltaX=(XMax-XMin)/DeltaX
    DeltaY=NumYDivisions%
    DeltaY=(YMax-YMin)/DeltaY
    X=XMin
    NumXPrimes%=0
    XDivisionNum%=1
    WHILE XDivisionNum% <= NumXDivisions%
      Y=YMin
      YDivisionNum%=1
      WHILE YDivisionNum% <= NumYDivisions%
        Z=F:(X,Y)
        NumXPrimes%=NumXPrimes%+1
        XDivisionIndex%(NumXPrimes%)=XDivisionNum%
        YDivisionIndex%(NumXPrimes%)=YDivisionNum%
        XRotated=X*CosRotation+Y*SinRotation
        YPrime(NumXPrimes%)=-X*SinRotation+Y*CosRotation
        XPrime(NumXPrimes%)=XRotated*CosTilt+Z*SinTilt
        ZPrime(NumXPrimes%)=-XRotated*SinTilt+Z*CosTilt
        IF XPrime(NumXPrimes%) > XPrimeMax
          XPrimeMax=XPrime(NumXPrimes%)
        ENDIF
        IF YPrime(NumXPrimes%) < YPrimeMin
          YPrimeMin=YPrime(NumXPrimes%)
        ENDIF
        IF YPrime(NumXPrimes%) > YPrimeMax
          YPrimeMax=YPrime(NumXPrimes%)
        ENDIF
        IF ZPrime(NumXPrimes%) < ZPrimeMin
          ZPrimeMin=ZPrime(NumXPrimes%)
        ENDIF
        IF ZPrime(NumXPrimes%) > ZPrimeMax
          ZPrimeMax=ZPrime(NumXPrimes%)
        ENDIF
        Y=Y+DeltaY
        YDivisionNum%=YDivisionNum%+1
      ENDWH
      X=X+DeltaX
      XDivisionNum%=XDivisionNum%+1
    ENDWH
    PRINT "Shading..."
    LightX=1.0
    LightY=-1.0
    LightZ=1.0
    Magnitude=SQR(LightX*LightX+LightY*LightY+LightZ*LightZ)
    LightX=LightX/Magnitude
    LightY=LightY/Magnitude
    LightZ=LightZ/Magnitude
    ColorMin%=255
    ColorMax%=0
    NumXPrimes%=0
    XDivisionNum%=1
    WHILE XDivisionNum% <= NumXDivisions%
      YDivisionNum%=1
      WHILE YDivisionNum% <= NumYDivisions%
        NumXPrimes%=NumXPrimes%+1
        VertexX(1)=XPrime(NumXPrimes%)
        VertexY(1)=YPrime(NumXPrimes%)
        VertexZ(1)=ZPrime(NumXPrimes%)
        IF XDivisionNum% < NumXDivisions%
          IF YDivisionNum% < NumYDivisions%
            NumXPrimes%=NumXPrimes%+NumYDivisions%
            VertexX(2)=XPrime(NumXPrimes%)
            VertexY(2)=YPrime(NumXPrimes%)
            VertexZ(2)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%+1
            VertexX(3)=XPrime(NumXPrimes%)
            VertexY(3)=YPrime(NumXPrimes%)
            VertexZ(3)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%-NumYDivisions%
            VertexX(4)=XPrime(NumXPrimes%)
            VertexY(4)=YPrime(NumXPrimes%)
            VertexZ(4)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%-1
          ELSE
            NumXPrimes%=NumXPrimes%-1
            VertexX(2)=XPrime(NumXPrimes%)
            VertexY(2)=YPrime(NumXPrimes%)
            VertexZ(2)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%+NumYDivisions%
            VertexX(3)=XPrime(NumXPrimes%)
            VertexY(3)=YPrime(NumXPrimes%)
            VertexZ(3)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%+1
            VertexX(4)=XPrime(NumXPrimes%)
            VertexY(4)=YPrime(NumXPrimes%)
            VertexZ(4)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%-NumYDivisions%
          ENDIF
        ELSE
          IF YDivisionNum% < NumYDivisions%
            NumXPrimes%=NumXPrimes%+1
            VertexX(2)=XPrime(NumXPrimes%)
            VertexY(2)=YPrime(NumXPrimes%)
            VertexZ(2)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%-NumYDivisions%
            VertexX(3)=XPrime(NumXPrimes%)
            VertexY(3)=YPrime(NumXPrimes%)
            VertexZ(3)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%-1
            VertexX(4)=XPrime(NumXPrimes%)
            VertexY(4)=YPrime(NumXPrimes%)
            VertexZ(4)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%+NumYDivisions%
          ELSE
            NumXPrimes%=NumXPrimes%-NumYDivisions%
            VertexX(2)=XPrime(NumXPrimes%)
            VertexY(2)=YPrime(NumXPrimes%)
            VertexZ(2)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%-1
            VertexX(3)=XPrime(NumXPrimes%)
            VertexY(3)=YPrime(NumXPrimes%)
            VertexZ(3)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%+NumYDivisions%
            VertexX(4)=XPrime(NumXPrimes%)
            VertexY(4)=YPrime(NumXPrimes%)
            VertexZ(4)=ZPrime(NumXPrimes%)
            NumXPrimes%=NumXPrimes%+1
          ENDIF
        ENDIF
        REM Compute the normal to a quadrilateral by averaging the
        REM normals to each vertex of the quadrilateral.
        NormalX=(VertexY(2)-VertexY(1))*(VertexZ(4)-VertexZ(1))
        NormalX=NormalX-(VertexY(4)-VertexY(1))*(VertexZ(2)-VertexZ(1))
        NormalX=NormalX+(VertexY(3)-VertexY(2))*(VertexZ(1)-VertexZ(2))
        NormalX=NormalX-(VertexY(1)-VertexY(2))*(VertexZ(3)-VertexZ(2))
        NormalX=NormalX+(VertexY(4)-VertexY(3))*(VertexZ(2)-VertexZ(3))
        NormalX=NormalX-(VertexY(2)-VertexY(3))*(VertexZ(4)-VertexZ(3))
        NormalX=NormalX+(VertexY(1)-VertexY(4))*(VertexZ(3)-VertexZ(4))
        NormalX=NormalX-(VertexY(3)-VertexY(4))*(VertexZ(1)-VertexZ(4))
        NormalY=(VertexX(4)-VertexX(1))*(VertexZ(2)-VertexZ(1))
        NormalY=NormalY-(VertexX(2)-VertexX(1))*(VertexZ(4)-VertexZ(1))
        NormalY=NormalY+(VertexX(1)-VertexX(2))*(VertexZ(3)-VertexZ(2))
        NormalY=NormalY-(VertexX(3)-VertexX(2))*(VertexZ(1)-VertexZ(2))
        NormalY=NormalY+(VertexX(2)-VertexX(3))*(VertexZ(4)-VertexZ(3))
        NormalY=NormalY-(VertexX(4)-VertexX(3))*(VertexZ(2)-VertexZ(3))
        NormalY=NormalY+(VertexX(3)-VertexX(4))*(VertexZ(1)-VertexZ(4))
        NormalY=NormalY-(VertexX(1)-VertexX(4))*(VertexZ(3)-VertexZ(4))
        NormalZ=(VertexX(2)-VertexX(1))*(VertexY(4)-VertexY(1))
        NormalZ=NormalZ-(VertexX(4)-VertexX(1))*(VertexY(2)-VertexY(1))
        NormalZ=NormalZ+(VertexX(3)-VertexX(2))*(VertexY(1)-VertexY(2))
        NormalZ=NormalZ-(VertexX(1)-VertexX(2))*(VertexY(3)-VertexY(2))
        NormalZ=NormalZ+(VertexX(4)-VertexX(3))*(VertexY(2)-VertexY(3))
        NormalZ=NormalZ-(VertexX(2)-VertexX(3))*(VertexY(4)-VertexY(3))
        NormalZ=NormalZ+(VertexX(1)-VertexX(4))*(VertexY(3)-VertexY(4))
        NormalZ=NormalZ-(VertexX(3)-VertexX(4))*(VertexY(1)-VertexY(4))
        Magnitude=SQR(NormalX*NormalX+NormalY*NormalY+NormalZ*NormalZ)
        IF Magnitude = 0.0
          ColorMin%=0
          Color%(NumXPrimes%)=0
        ELSE
          ColorNum%=INT((256.0/2.0)*(1.0+(LightX*NormalX+LightY*NormalY+LightZ*NormalZ)/Magnitude))
          IF ColorNum% > 255
            ColorNum%=255
          ENDIF
          Color%(NumXPrimes%)=ColorNum%
          IF ColorNum% < ColorMin%
            ColorMin%=ColorNum%
          ENDIF
          IF ColorNum% > ColorMax%
            ColorMax%=ColorNum%
          ENDIF
        ENDIF
        YDivisionNum%=YDivisionNum%+1
      ENDWH
      XDivisionNum%=XDivisionNum%+1
    ENDWH
    PRINT "Adjusting perspective..."
    IF YPrimeMax-YPrimeMin > ZPrimeMax-ZPrimeMin
      XEye=2.0*(YPrimeMax-YPrimeMin)+XPrimeMax
    ELSE
      XEye=2.0*(ZPrimeMax-ZPrimeMin)+XPrimeMax
    ENDIF
    IF ZPrimeMax <> ZPrimeMin
      YCenter=(YPrimeMax+YPrimeMin)/2.0
      ZCenter=(ZPrimeMax+ZPrimeMin)/2.0
      NumXPrimes%=0
      XDivisionNum%=1
      WHILE XDivisionNum% <= NumXDivisions%
        Y=YMin
        YDivisionNum%=1
        WHILE YDivisionNum% <= NumYDivisions%
          NumXPrimes%=NumXPrimes%+1
          X=XPrime(NumXPrimes%)
          Y=YPrime(NumXPrimes%)
          Z=ZPrime(NumXPrimes%)
          DeltaX=X-XEye
          DeltaY=Y-YCenter
          DeltaZ=Z-ZCenter
          XPrime(NumXPrimes%)=SQR(DeltaX*DeltaX+DeltaY*DeltaY+DeltaZ*DeltaZ)
          YPrime(NumXPrimes%)=YCenter+(Y-YCenter)*(XEye-XPrimeMax)/(XEye-X)
          ZPrime(NumXPrimes%)=ZCenter+(Z-ZCenter)*(XEye-XPrimeMax)/(XEye-X)
          YDivisionNum%=YDivisionNum%+1
        ENDWH
        XDivisionNum%=XDivisionNum%+1
      ENDWH
    ENDIF
    PRINT "Sorting points..."
    SortLeft%=NumXPrimes%/2
    SortLeft%=SortLeft%+1
    SortRight%=NumXPrimes%
    SortT1=XPrime(1)
    SortT2%=XDivisionIndex%(1)
    SortT3%=YDivisionIndex%(1)
    WHILE SortRight% > 1
      IF SortLeft% > 1
        SortLeft%=SortLeft%-1
        SortT1=XPrime(SortLeft%)
        SortT2%=XDivisionIndex%(SortLeft%)
        SortT3%=YDivisionIndex%(SortLeft%)
      ELSE
        SortT1=XPrime(SortRight%)
        SortT2%=XDivisionIndex%(SortRight%)
        SortT3%=YDivisionIndex%(SortRight%)
        XPrime(SortRight%)=XPrime(1)
        XDivisionIndex%(SortRight%)=XDivisionIndex%(1)
        YDivisionIndex%(SortRight%)=YDivisionIndex%(1)
        SortRight%=SortRight%-1
      ENDIF
      IF SortRight% > 1
        KeyIndex2%=SortLeft%
        Finished%=0
        WHILE Finished% = 0
          KeyIndex1%=KeyIndex2%
          KeyIndex2%=2*KeyIndex2%
          IF KeyIndex2% > SortRight%
            Finished%=-1
          ELSE
            IF KeyIndex2% <> SortRight%
              IF XPrime(KeyIndex2%) > XPrime(KeyIndex2%+1)
                KeyIndex2%=KeyIndex2%+1
              ENDIF
            ENDIF
            IF SortT1 <= XPrime(KeyIndex2%)
              Finished%=-1
            ELSE
              XPrime(KeyIndex1%)=XPrime(KeyIndex2%)
              XDivisionIndex%(KeyIndex1%)=XDivisionIndex%(KeyIndex2%)
              YDivisionIndex%(KeyIndex1%)=YDivisionIndex%(KeyIndex2%)
            ENDIF
          ENDIF
        ENDWH
        XPrime(KeyIndex1%)=SortT1
        XDivisionIndex%(KeyIndex1%)=SortT2%
        YDivisionIndex%(KeyIndex1%)=SortT3%
      ENDIF
    ENDWH
    XPrime(1)=SortT1
    XDivisionIndex%(1)=SortT2%
    YDivisionIndex%(1)=SortT3%
    idPlot%=GCREATE(0,0,640,240,1,2)
    REM AspectRatio=1.0/(4.0*(240.0/640.0)/3.0)
    AspectRatio=1.0/(54.0*(240.0/640.0)/21.0)
    YOutMax=MaxYOut%
    ZOutMax=MaxZOut%
    IF AspectRatio*ZOutMax*(YPrimeMax-YPrimeMin) <= YOutMax*(ZPrimeMax-ZPrimeMin)
      IF AspectRatio*ZOutMax*(YPrimeMax-YPrimeMin) >= YOutMax*(ZPrimeMax-ZPrimeMin)
        PixelsPerUnit=1.0
        YOffset=YOutMax/2.0
        ZOffset=-ZOutMax/2.0
      ELSE
        PixelsPerUnit=ZOutMax/(ZPrimeMax-ZPrimeMin)
        YOffset=(YOutMax-AspectRatio*PixelsPerUnit*(YPrimeMax-YPrimeMin))/2.0
        ZOffset=0.0
      ENDIF
    ELSE
      PixelsPerUnit=YOutMax/(AspectRatio*(YPrimeMax-YPrimeMin))
      YOffset=0.0
      ZOffset=-(ZOutMax-PixelsPerUnit*(ZPrimeMax-ZPrimeMin))/2.0
    ENDIF
    Bias=1.0
    LightestGray%=240
    DarkestGray%=19
    XPrimeNum%=1
    WHILE XPrimeNum% <= NumXPrimes%
      XDivisionNum%=XDivisionIndex%(XPrimeNum%)
      IF XDivisionNum% <> NumXDivisions%
        YDivisionNum%=YDivisionIndex%(XPrimeNum%)
        IF YDivisionNum% <> NumYDivisions%
          DivIndex%=DivisionIndex%:(XDivisionNum%,YDivisionNum%,NumXDivisions%,NumYDivisions%)
          BoxX%(1)=INT(YOffset+PixelsPerUnit*AspectRatio*(YPrime(DivIndex%)-YPrimeMin))
          BoxY%(1)=INT(ZOffset+ZOutMax-PixelsPerUnit*(ZPrime(DivIndex%)-ZPrimeMin))
          BoxX%(2)=INT(YOffset+PixelsPerUnit*AspectRatio*(YPrime(DivIndex%+NumYDivisions%)-YPrimeMin))
          BoxY%(2)=INT(ZOffset+ZOutMax-PixelsPerUnit*(ZPrime(DivIndex%+NumYDivisions%)-ZPrimeMin))
          BoxX%(3)=INT(YOffset+PixelsPerUnit*AspectRatio*(YPrime(DivIndex%+NumYDivisions%+1)-YPrimeMin))
          BoxY%(3)=INT(ZOffset+ZOutMax-PixelsPerUnit*(ZPrime(DivIndex%+NumYDivisions%+1)-ZPrimeMin))
          BoxX%(4)=INT(YOffset+PixelsPerUnit*AspectRatio*(YPrime(DivIndex%+1)-YPrimeMin))
          BoxY%(4)=INT(ZOffset+ZOutMax-PixelsPerUnit*(ZPrime(DivIndex%+1)-ZPrimeMin))
          BoxYMin%=BoxY%(1)
          BoxYMax%=BoxYMin%
          BoxNum1%=2
          WHILE BoxNum1% <= 4
            IF BoxY%(BoxNum1%) < BoxYMin%
              BoxYMin%=BoxY%(BoxNum1%)
            ENDIF
            IF BoxY%(BoxNum1%) > BoxYMax%
              BoxYMax%=BoxY%(BoxNum1%)
            ENDIF
            BoxNum1%=BoxNum1%+1
          ENDWH
          ColorNum%=Color%(DivIndex%)
          Fraction=FLT(ColorNum%-ColorMin%)/FLT(ColorMax%-ColorMin%)
          IF Fraction > 0.0
            Fraction=EXP(Bias*LN(Fraction))
            ColorNum%=INT(FLT(LightestGray%-DarkestGray%)*Fraction+FLT(DarkestGray%))
          ELSE
            ColorNum%=DarkestGray%
          ENDIF
          GCOLOR ColorNum%,ColorNum%,ColorNum%
          BoxY1%=BoxYMin%
          WHILE BoxY1% <= BoxYMax%
            InterceptCountMod2%=0
            BoxNum2%=2
            BoxNum1%=1
            WHILE BoxNum1% <= 4
              IF BoxY%(BoxNum1%) >= BoxY1%
                IF BoxY1% > BoxY%(BoxNum2%)
                  BoxDeltaY=BoxY%(BoxNum2%)-BoxY%(BoxNum1%)
                  BoxDeltaX=BoxX%(BoxNum2%)-BoxX%(BoxNum1%)
                  BoxYOffset=BoxY1%-BoxY%(BoxNum1%)
                  BoxXIntercept=BoxX%(BoxNum1%)
                  BoxX1%=INT(BoxDeltaX*BoxYOffset/BoxDeltaY+BoxXIntercept)
                  DrawLine:
                  InterceptCountMod2%=1-InterceptCountMod2%
                ENDIF
              ELSE
                IF BoxY1% <= BoxY%(BoxNum2%)
                  BoxDeltaY=BoxY%(BoxNum2%)-BoxY%(BoxNum1%)
                  BoxDeltaX=BoxX%(BoxNum2%)-BoxX%(BoxNum1%)
                  BoxYOffset=BoxY1%-BoxY%(BoxNum1%)
                  BoxXIntercept=BoxX%(BoxNum1%)
                  BoxX1%=INT(BoxDeltaX*BoxYOffset/BoxDeltaY+BoxXIntercept)
                  DrawLine:
                  InterceptCountMod2%=1-InterceptCountMod2%
                ENDIF
              ENDIF
              BoxNum2%=BoxNum2%+1
              IF BoxNum2% > 4
                BoxNum2%=1
              ENDIF
              BoxNum1%=BoxNum1%+1
            ENDWH
            BoxY1%=BoxY1%+1
          ENDWH
        ENDIF
      ENDIF
      XPrimeNum%=XPrimeNum%+1
    ENDWH
    Response$=GET$
    GCLOSE idPlot%
    CLS
    PRINT "                             Three Dimensional Plot"
    PRINT
    PRINT
    PRINT
    PRINT "Again (y or n)? ";
    Response$=GET$
    PRINT
  UNTIL Response$ <> "Y" AND Response$ <> "y"
ENDP

PROC DrawLine:
  EXTERNAL BoxX1%
  EXTERNAL BoxX2%
  EXTERNAL BoxY1%
  EXTERNAL InterceptCountMod2%
  IF InterceptCountMod2% = 0
    BoxX2%=BoxX1%
  ELSE
    GAT BoxX1%,BoxY1%
    GLINETO BoxX2%,BoxY1%
    GLINEBY 0,0
  ENDIF
ENDP

PROC DivisionIndex%:(XDivisionNum%,YDivisionNum%,NumXDivisions%,NumYDivisions%)
  RETURN NumYDivisions%*(XDivisionNum%-1)+YDivisionNum%
ENDP

PROC F:(X,Y)
  REM     This function plots a circular wave in three dimensions.
  REM Try plotting -1 <= x <= 1 and -1 <= y <= 1 with 30 divisions
  REM along each axis, a rotation of 30 degrees, and a tilt of 20
  REM degrees.
  LOCAL T1
  LOCAL T2
  T1=X*X+Y*Y
  T2=COS(7.0*SQR(T1))
  RETURN 2*T2*T2/(1.0+30.0*T1)
ENDP