Example showing drawing using GINOSURF

!      Program SURFEX7
! **********************************************
! Copyright Bradly Associates Limited
!        GINOSURF  Version 6.0
! **********************************************
!
!  projected contour map and 4d triangulated surface
!  on 3d axes using computed data
!
! *************************************************
!
use surf_f90
!
!  declare variables
!
parameter (MAXPOINTS=1000, NCONT=24)
!
type (GWORKSPACE)    work
!
type (GPOINT4) :: pts(MAXPOINTS),bnd(100)
type (GPOINT)  :: deriv(1),deriv2(1)
type (GDIM)       paper
type (GDEVSTATE)  device_state
type (GLIMIT)  :: picture2D = GLIMIT(0.0,1000.0,0.0,750.0)
type (GLIMIT3) :: picture3D = GLIMIT3(0.0,1000.0,0.0,750.0,-5000.0,20000.0)
type (GLIMIT)  :: window1 = GLIMIT(200.0,900.0,100.0,700.0)
type (GLIMIT)  :: frame1 = GLIMIT(10.0,990.0,50.0,740.0)
type (GLIMIT)  :: frame2 = GLIMIT(15.0,985.0,55.0,735.0)
type (GLIMIT)  :: data_limit = GLIMIT(200.0,500.0,100.0,340.0)
type (GLIMIT)  :: keybox = GLIMIT(80.0,120.0,100.0,400.0)
type (GLIMIT)  :: viewport = GLIMIT(0.0,1.0,0.0,1.0)
type (GLIMIT3)    data_3D
type (GFILL)      index(NCONT+1)
type (GSURF3D) :: state = GSURF3D(GON,GOFF)
type (GPRJAXI) :: style = GPRJAXI(GAXISANDGRID,GAXISANDGRID,GAXISANDGRID,GAXISANDGRID,GAXISANDGRID,GAXISANDGRID)
type (GPRJAXI) :: supress = GPRJAXI(GNONE,GNONE,GNONE,GNONE,GNONE,GNONE)

      real zlev(NCONT+1)
!
!  Initialise variables
!
      zmin = 999999.9
      zmax = -zmin
      np = 0
      nb = 0
      nbmax = 72
!
!  Generate coordinates, both z values and hole at the pole
!
      pi180 = 3.14159265359/180.0
      do irad = 2, 10
         rad = 10.0*sin(irad*pi180*9.0)
         do iang = 5, 180, 5
            ang = iang * pi180
            np = np + 1
            pts(np)%x = rad*cos(ang)
            pts(np)%y = rad*sin(ang)
            if (abs(pts(np)%x).lt.1.0e-5) pts(np)%x = 0.0
            if (abs(pts(np)%y).lt.1.0e-5) pts(np)%y = 0.0
            pts(np)%z = sqrt(abs(100.0- pts(np)%x**2 - pts(np)%y**2))
            pts(np)%z2 = sqrt(abs(80.0-(pts(np)%x)**2-5*(pts(np)%y+5)**2))
            if (pts(np)%z2 .gt. zmax) zmax = pts(np)%z2
            if (pts(np)%z2 .lt. zmin) zmin = pts(np)%z2
            if (irad .eq. 2) then
               nb = nb + 1
               bnd(nb)%x = pts(np)%x
               bnd(nb)%y = pts(np)%y
               bnd(nb)%z = pts(np)%z
               bnd(nb)%z2 = pts(np)%z2
            end if
            if(pts(np)%y.ne.0.0) then
               np = np + 1
               pts(np)%x = pts(np-1)%x
               pts(np)%y = -pts(np-1)%y
               pts(np)%z = pts(np-1)%z
               pts(np)%z2 = sqrt(abs(80.0-(pts(np)%x)**2-5*(pts(np)%y+5)**2))
               if (pts(np)%z2 .gt. zmax) zmax = pts(np)%z2
               if (pts(np)%z2 .lt. zmin) zmin = pts(np)%z2
            end if
            if (irad .eq. 2) then
               nbt = nbmax - nb + 1
               bnd(nbt)%x = pts(np)%x
               bnd(nbt)%y = pts(np)%y
               bnd(nbt)%z = pts(np)%z
               bnd(nbt)%z2 = pts(np)%z2
            end if
         end do
      end do
!
!
!  Initialize gino
!
      call gOpenGino
      call xxxxx
      call gEnqDrawingLimits(paper,ip)
!
!  Check for 3D device
!
      call gEnqDeviceState(device_state)
      viewport%xmax=paper%xpap
      viewport%ymax=paper%ypap
      if(device_state%dddim .eq. 2) then
         call gSetViewport2D(picture2D,viewport)
      else
         call gSetViewport3D(picture3D,viewport)
      end if
      call gNewDrawing
      call gSetCharSize(8.5,8.5)
      call gEnqColourInfo(ndc,ndt)
      call gEnqLineColour(icol)
!
!  Define colour table
!
      call gDefineRGB(2,0.4,0.0,0.0)
      call gDefineRGB(3,1.0,0.0,0.0)
      call gDefineRGB(4,1.0,0.3,0.3)
      call gDefineRGB(5,0.9,0.75,0.1)
      call gDefineRGB(6,1.0,0.9,0.1)
      call gDefineRGB(7,0.9,1.0,0.3)
      call gDefineRGB(8,0.1,1.0,0.4)
      call gDefineRGB(9,0.0,0.8,0.5)
      call gDefineRGB(10,0.0,0.95,0.85)
      call gDefineRGB(11,0.0,0.2,1.0)
      call gDefineRGB(12,0.4,0.4,0.8)
      call gDefineRGB(13,0.7,0.5,0.9)
      call gDefineRGB(14,0.5,0.1,0.9)
      call gDefineRGB(15,0.3,0.3,0.3)
!
!  Initialise surface drawing
!
      call gsOpenSurf
      call gsSetSurf3DState(state)
      call gsSetContourMapFrame(window1)
      call gsSetSurfaceFrame(window1)
!
!  Generate 4d triangular network from random points
!
      call gsGenerateRandomNetwork4D(MAXPOINTS,np,pts,work)
!
!  Extract central region
!
      call gsAddNetworkRegion4D(GXYZONLY,nbmax,bnd,deriv,deriv2,GEXTERIOR,work)
!
!  Define contour levels
!
      zmin = 0.0
      zstep = 1.5
      zlev(1) = zmin
      do i = 2, NCONT+1
         zlev(i) = zmin + zstep*(i-1)
      end do
      call gsSetContourLevels(NCONT-1,zlev)
      call gsSetContourDrawingSwitch(GON)
!
!  Define contour key box attrbutes
!
      key = GON
      call gsSetContourKeyPos(keybox,gsKeyend=GROUNDED)
!
!  Define fill styles
!
      do i = 1, NCONT+1
         it = i + 1
         if(it.gt.15) it = it - 14
         index(i)%fill=GSOLID
         index(i)%line=it
      end do
!
!  Define perspective view
!
      data_limit%xmin = -10.0
      data_limit%ymin = -10.0
      data_limit%xmax = 10.0
      data_limit%ymax = 10.0
      data_3D%xmin=data_limit%xmin
      data_3D%xmax=data_limit%xmax
      data_3D%ymin=data_limit%ymin
      data_3D%ymax=data_limit%ymax
      data_3D%zmin = -20.0
      data_3D%zmax = 10.0
      radius = 500.0
      theta = 25.0
      phi = 30.0
      call gsSetSurfaceBaseType(GBOTHSURFACES)
      call gsSetSurfaceAxesLabelMode(GADJACENT)
      call gsSetSurfaceAxesStyle(style,supress)
      call gsDefinePerspProjection(data_3D, radius, theta, phi, work)
      call gsDrawSurfaceAxes(work)
!
!  Specify to contour from z2 and set display height
!
      call gsSelectNetworkContourDataSet(GSECOND)
      call gsSetContourMapHeight(GON,-20.0)
      call gsFillTriangulatedContourMap(data_limit,NCONT,GSTRAIGHT,index,key,work)
!
!  Add 4d perspective surface
!
      call gsFillTriangulatedContourSurf4D(data_limit, radius, theta, phi, NCONT, index, key, work)
!
!  Add 3d axes and triangulation network
!
      call gsDrawNetwork(GNONE,0,work)
!
!  Close down GINOSURF
!
      call gsFreeWorkspace(work)
      call gsCloseSurf
!
!  Draw frame and title
!
      call gSetWindowMode(GON2D)
      call gSetLineColour(GWHITE)
      call gSetCharSize(15.0,15.0)
      call gFillRect(GHOLLOW,GCURRENT,frame1)
      call gFillRect(GHOLLOW,GCURRENT,frame2)
      call gMoveTo2D(305.0,20.0)
      call gDisplayStr('GINOSURF Example Program 7')
      call gSetStrJustify(GCENTRE)
      call gMoveTo2D(500.0,700.0)
      call gDisplayStr('4D Triangulated Surface and Projected Contour')
!
!  Close down GINO
!
      call gSuspendDevice
      call gCloseGino
      stop
      end