' ******************************************************************
' Contour plot of a two-dimensional function
' ******************************************************************

#include "fbmath.bi"
#include "fbgfx.bi"

' ******************************************************************

' ******************************************************************
' Define here the function to be plotted
' ******************************************************************

CONST FuncName = "SIN(SQR(X^2 + Y^2)) + 0.5 / SQR((X + 3.05)^2 + Y^2)"

FUNCTION Func(X AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Function to be plotted

  CONST C = 3.05

  DIM AS DOUBLE X2, Y2, Xc, Xc2

  X2 = X * X
  Y2 = Y * Y
  Xc = X + C
  Xc2 = Xc * Xc

  Func = SIN(SQR(X2 + Y2)) + 0.5 / SQR(Xc2 + Y2)
END FUNCTION

' ******************************************************************
' Main program
' ******************************************************************

' Initialize screen (Screen 12 with white on black here)

IF NOT InitScreen(12, 15, 0) THEN STOP

' Set number of contours (colours) according to the screen resolution

CONST Nc = 15

' Define colours

DIM AS INTEGER Col(0 TO Nc - 1), I

FOR I = 0 TO Nc - 1
  Col(I) = Nc - I
NEXT I

' Set scale

CONST Xmin = - TwoPi, Xmax = TwoPi, Xstep = Pi
CONST Ymin = - TwoPi, Ymax = TwoPi, Ystep = Pi

' Set graphic area (here it will look approximately square)

InitGraph 24, 76, 15, 85

' Set number of grid points

CONST NpX = 60
CONST NpY = 60

' Dimension variables

DIM AS INTEGER J                           ' Loop variable
DIM AS DOUBLE  Dx, Dy                      ' Increments of X and Y
DIM AS DOUBLE  Dmin, Dmax                  ' Min. and max. function values
DIM AS DOUBLE  H                           ' Increment for levels

DIM AS DOUBLE  Xu(0 TO NpX), Yu(0 TO NpY)  ' Point coordinates in user units
DIM AS INTEGER Xp(0 TO NpX), Yp(0 TO NpY)  ' Point coordinates in pixels
DIM AS DOUBLE  D(0 TO NpX, 0 TO NpY)       ' Function values
DIM AS DOUBLE  Z(0 TO Nc - 1)              ' Contour array

' Set scales and plot axes

SetOxScale LinScale, Xmin, Xmax, Xstep
SetOyScale LinScale, Ymin, Ymax, Ystep

PlotOxAxis Ymin, 2
PlotOyAxis Xmin, 2

PlotGrid BothGrid

WriteTitles FuncName, "X", "Y"

' Generate grid points

Dx = (Xmax - Xmin) / NpX
Dy = (Ymax - Ymin) / NpY

Xu(0) = Xmin
FOR I = 1 TO NpX
  Xu(I) = Xu(I - 1) + Dx
NEXT I

Yu(0) = Ymin
FOR J = 1 TO NpY
  Yu(J) = Yu(J - 1) + Dy
NEXT J

FOR I = 0 TO NpX
  Xp(I) = Xpixel(Xu(I))
NEXT I

FOR J = 0 TO NpY
  Yp(J) = Ypixel(Yu(J))
NEXT J

' Compute function values

Dmin = Func(Xmin, Ymin)
Dmax = Dmin

FOR I = 0 TO NpX
  FOR J = 0 TO NpY
    D(I,J) = Func(Xu(I), Yu(J))
    IF D(I,J) < Dmin THEN
      Dmin = D(I,J)
    ELSEIF D(I,J) > Dmax THEN
      Dmax = D(I,J)
    END IF
  NEXT J
NEXT I

' Define levels

Z(0) = Dmin

H = (Dmax - Dmin) / (Nc - 1)

FOR I = 1 TO Nc - 1
  Z(I) = Z(I - 1) + H
NEXT I

' Plot contour

ConRec Xp(), Yp(), Z(), D(), Col()

SLEEP
END
