c     *******************************************
c     * Terrain Geometry Generation Program
c     *
c     * Kent F. Martin
c     * 15 JAN 90
c     * Animated Technologies
c     * (213) 675 - 0770
c     * 
c     * This program reads Defense Terrain Mapping data and
c     * converts is to a geometry file format which is compatible
c     * with WaveFront Technology software (Model, Preview, Image)
c     *
c     * This version offers the choice of high, normal, or low
c     * resolution output. Normal resolution output has polygons
c     * with vertices located at the vertices in the database,
c     * high resolution adds one vertex between the original
c     * ones, resulting in 4 times more polygons, low resolution
c     * skips some data values depending on user choice as to how 
c     * much data to skip between saves, this results in considerably
c     * less polygons
c     *
c     ******************************************

c.....Variable description list

c.....arrays:
c       lomin : minimum longitude of database
c       lomax : maximum longitude of database
c       lamin : minimum latitude of database
c       lamax : maximum latitude of database
c
c       lolo : minimum longitude of geometry region
c       uplo : maximum longitude of geometry region
c       lola : minimum latitude of geometry region
c       upla : maximum latitude of geometry region
c
c       a1 : storage array #1, holds elevation data for processing
c       a2 : storage array #2, holds elevation data for processing
c            these are the only two storage arrays, making possible
c            the processing of a large database without the need
c            for large memory for program initialization.
c
c       losave : storage array for data when longitude values are being
c                skipped in lo resolution, serves same purpose as a1 array
c.....single variables
c
c     v1,v2,v3,...,v9 : temporary vertex numbers for polygon creation
c
c     xint : x coordinate midway between lo and (lo-1)
c     yint : y coordinate midway between la and (la-1)
c     int  : elevation at a point utilizing xint or yint, but not both
c     int1 : elevation at a point when two intermediate points exist
c     int2 : elevation at a point when two intermediate points exist
c
c
c     lotot : counter of number of longitude values in database
c     latot : number of latitude values in database
c     lost : counter where processing of longitude data begins
c     loend: counter where processing of longitude data ends
c     last : counter where processing of latitude data begins
c     laend: counter where processing of latitude data ends

c     ic   : integer seed value for intermediate point altitude
c            generation
c     z    : intermediate point altitude, result of using then
c            Gaussian distribution random number generator
c     ave  : average of elevations of points surrounding 
c            intermediate point
c     skip : skip factor for low resolution geometry
c
c     lo : longitude counter
c     la : latitude counter
c
c.....file names
c     in : file containing elevation data
c     out: output file, this will be in a form which "model"
c          can read, suggest user supply file.obj for name.

      character*80 in, out
      real a1(5000), a2(5000)
      real int, int1, int2
      integer lolo(3), uplo(3), lola(3), upla(3)
      integer lomax(3), lomin(3), lamax(3), lamin(3)
      integer v1,v2,v3,v4,v5,v6,v7,v8,v9
      integer skip
      integer losave
      integer choice

      print*,'Welcome to the Terrain Generator'
      print*,' '
      print*,'Enter data file name to be processed'
      print*,' '
      read(*,*)in
      print*,'Enter output object file name'
      print*,' '
      read(*,*)out
      print*,' '
      print*,'Select desired level of terrain resolution:'
      print*,'------'
      print*,'  1) high'
      print*,'  2) normal'
      print*,'  3) low'
      print*,' '
      read(*,*)level

c.....if low level, ask for enlargement factor

      if ( level .eq. 3 ) then
        print*,' '
        print*,'Enter enlargement factor for geometry (integer)'
        read(*,*)skip
      endif

c.....read lower and upper limits for geometry boundary

      print*,'This program is intended to be used for elevation'
      print*,'data recorded at 3 second intervals of longitude and'
      print*,'latitude. If the value for seconds(ss) is not an'
      print*,'integer multiple of 3, the input will be requested'
      print*,'again.'

c.....enter database limits of longitude and latitude
 6    print*,' '
      print*,'Enter database lower longitude limit ( dd(d) mm ss )'
      read(*,*)(lomin(i),i=1,3)
      if(mod(lomin(3),3) .ne. 0 ) goto 6

 7    print*,' '
      print*,'Enter database upper longitude limit ( dd(d) mm ss )'
      read(*,*)(lomax(i),i=1,3)
      if(mod(lomax(3),3) .ne. 0 ) goto 7

 8    print*,' '
      print*,'Enter database lower latitude limit ( dd mm ss )'
      read(*,*)(lamin(i),i=1,3)
      if(mod(lamin(3),3) .ne. 0 ) goto 8

 9    print*,' '
      print*,'Enter database upper latitude limit ( dd mm ss )'
      read(*,*)(lamax(i),i=1,3)
      if(mod(lamax(3),3) .ne. 0 ) goto 9

c.....specify area of database to be processed

      print*,' '
 10   print*,'Enter lower longitude limit ( dd(d) mm ss )'
      read(*,*)(lolo(i),i=1,3)
      if( mod(lolo(3),3) .ne. 0 )goto 10

      print*,' '
 11   print*,'Enter upper longitude limit ( dd(d) mm ss )'
      read(*,*)(uplo(i),i=1,3)
      if( mod(uplo(3),3) .ne. 0 )goto 11

      print*,' '
 12   print*,'Enter lower latitude limit ( dd mm ss )'
      read(*,*)(lola(i),i=1,3)
      if( mod(lola(3),3) .ne. 0)goto 12

      print*,' '
 13   print*,'Enter upper latitude limit ( dd mm ss )'
      read(*,*)(upla(i),i=1,3)
      if( mod(upla(3),3) .ne. 0)goto13

c....., form='unformated'open input and output files

      open(unit=20, file=in, status='old',form='unformatted',
     & access='direct',recl=9604)
      open(unit=21, file=out, status='unknown')

c.....establish do loop control counters by using geodiff subroutine
c.....minimum and maximum longitude and latitude should be read in from
c.....data file header, if not it must be requested interactively
c.....they will be hard coded for this application

      call geodiff(lolo, lomin, lost)
      call geodiff(uplo, lomin, loend)
      call geodiff(lola, lamin, last)
      call geodiff(upla, lamin, laend)
      call geodiff(lomax, lomin, lotot)
      call geodiff(lamax, lamin, latot)
c.....print control counters for debugging
c     print*,'lotot =',lotot
c     print*,'latot =',latot
c     print*,'lost =',lost
c     print*,'loend =',loend
c     print*,'last =',last
c     print*,'laend =',laend

c.....compute number of vertices and polygons
c.....notify user of these
c.....give option to reenter limits, continue, or quit.

      if(level.eq.1)then
        numv=((loend-lost+1)*2-1)*((laend-last+1)*2-1)
        nump=(loend-lost)*(laend-last)*8
      elseif(level.eq.2)then
        numv=(loend-lost+1)*(laend-last+1)
        nump=(loend-lost)*(laend-last)*2
      elseif(level.eq.3)then
        if(mod(loend-lost,skip) .eq. 0 )then
          loadd = 1
        else
          loadd = 2
        endif
        if(mod(laend-last,skip) .eq. 0 )then
          laadd = 1
        else
          laadd = 2
        endif
        numv=((loend-lost)/skip+loadd)*((laend-last)/skip+laadd)
        nump=((loend-lost)/skip+loadd-1)*
     &       ((laend-last) /skip+laadd-1) * 2
      endif

      print*,' '
      print*,'Current geography limits and resolution type'
      print*,'        will yield:'
      print*,' ',numv,' vertices and'
      print*,' ',nump,' polygons.'
      print*,' '
 80   print*,'Do you wish to:'
      print*,'  1) continue'
      print*,'  2) enter new geography limits'
      print*,'  3) quit?'
      read(*,*)choice

      if(choice .eq. 1)then
        goto 90
      elseif(choice .eq. 2)then
        goto 6
      elseif(choice .eq. 3)then
        goto 1000
      else
        goto 80
      endif

c.....branch to correct area depending on resolution level

c.....high resolution area

 90   continue
      if ( level .eq. 1 ) then

c.....do only enough lines to cover longitude range

        do 100 lo = 1, loend

          if ( lo .lt. lost ) then
c..... a dummy line if longitude out of range
            read(20,rec=lo)dummy
c.....            format(a1)
          elseif ( lo .eq. lost ) then
c.....read first line of data and write vertex information
            read(20,rec=lo)(a1(la),la=1,latot)
            do 110 la = last, laend
              if ( la .eq. last ) then
c.....write only one point if latitude is minimum value
                x = 60. * ( lo - 1)
                y = 60. * ( la - 1)
                write(21,*)'v ',x,y,a1(la)
              else
c.....write two point if latitude is greater than minimum
                x = 60. * ( lo -1)
                y = 60. * ( la -1)
                yint = y - 30.
c.....intermediate grid point #1
                ic = lo * ( la - 1 )
                z = grn(ic)
                ave = (a1(la) + a1(la-1)) / 2.
                deltaz = abs( a1(la) - a1(la-1) )
                int = ave + deltaz / 10. * z
                write(21,*)'v ',x,yint,int
                write(21,*)'v ',x,y,a1(la)
              endif
 110        continue
          else
c.....if longitude is not starting value, ad a line of data
c.....and store in second array
           read(20,rec=lo)(a2(la),la=1,latot)
           do 200 la = last, laend
             if ( la .eq. last ) then
c.....write two vertices for minimum latitude
               xint = 60.0 * ( lo - 1) - 30.
               y = 60. * (la -1)
c.....intermediate grid point #2
               ic = (lo - 1) * ( la - 1)
               z = grn(ic)
               ave = ( a1(la) + a2(la) ) / 2.
               deltaz = abs ( a1(la) - a2(la) )
               int = ave + deltaz / 10. * z
               write(21,*)'v ',xint,y,int
             else
c.....determine and write two more intermediate points
c.....intermediate grid point #3
               ic = (lo-1)*(la-1)+(lo*la)
               z = grn(ic)
               ave =(a1(la)+a2(la)+a1(la-1)+a2(la-1))/4.
               zmax = max(a1(la),a2(la),a1(la-1),a2(la-1))
               zmin = min(a1(la),a2(la),a1(la-1),a2(la-1))
               deltaz = abs(zmax-zmin)
               int1 = ave + deltaz / 10. * z
c.....intermediate grid point #4
               ic = (lo-1)*la
               z = grn(ic)
               ave = (a1(la) + a2(la)) / 2.
               deltaz = abs( a1(la) - a2(la) )
               int2 = ave + deltaz / 10. * z
               xint = 60. * ( lo -1) - 30.
               yint = 60. * ( la -1) - 30.
               y = 60. * (la -1)
               write(21,*)'v ',xint,yint,int1
               write(21,*)'v ',xint,y,int2
             endif
 200       continue

c.....without increasin longitude counter, go through latitude
c.....values generating the right hand side of the geometry 
c.....generated to this point.
           do 300 la = last, laend
             if ( la .eq. last ) then
               x = 60. * ( lo -1)
               y = 60. * ( la -1)
               write(21,*)'v ',x,y,a2(la)
             else
               x = 60. * ( lo -1)
               y = 60. * ( la -1)
               yint = y - 30.
c.....intermediate point #5
               ic = lo * la
               z = grn(ic)
               ave = ( a2(la) + a2(la-1) ) / 2.
               deltaz = abs ( a2(la) - a2(la-1) )
               int = ave + deltaz / 10. * z
               write(21,*)'v ',x,yint,int
               write(21,*)'v ',x,y,a2(la)
             endif
 300       continue

c.....begin the creation of polygons
c.....start with a longitude value one step larger than minimum

c.....polygons are built by analyzing four points in the
c.....original database which make a rectangle. this region
c.....is converted into 8 triangles by use of the intermediate
c.....points generated earlier. the 'direction' of the triangles
c.....is determined by finding the highest corner and working
c.....down from that point.
           do 400 la = last+1, laend
             e1 = a1(la-1)
             e2 = a2(la-1)
             e3 = a2(la)
             e4 = a1(la)
             print*,'lo =',lo
             print*,'la =',la
             print*,'lost =',lost
             print*,'last =',last
             print*,'laend =',laend
             v1=(lo-lost-1)*(((laend-last+1)*2-1)*2)+(la-last)*2 - 1
             v2 = v1 + 1
             v3 = v1 + 2
             v4 = v1 + (((laend-last+1)*2)-1)
             v5 = v4 + 1
             v6 = v4 + 2
             v7 = v4 + (((laend-last+1)*2)-1)
             v8 = v7 + 1
             v9 = v7 + 2
             if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
     &          e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
               write(21,*)'fo ',v1,v4,v5
               write(21,*)'fo ',v1,v5,v2
               write(21,*)'fo ',v2,v5,v6
               write(21,*)'fo ',v2,v6,v3
               write(21,*)'fo ',v4,v7,v8
               write(21,*)'fo ',v4,v8,v5
               write(21,*)'fo ',v5,v8,v9
               write(21,*)'fo ',v5,v9,v6
             else
               write(21,*)'fo ',v2,v1,v4
               write(21,*)'fo ',v2,v4,v5
               write(21,*)'fo ',v3,v2,v5
               write(21,*)'fo ',v3,v5,v6
               write(21,*)'fo ',v5,v4,v7
               write(21,*)'fo ',v5,v7,v8
               write(21,*)'fo ',v6,v5,v8
               write(21,*)'fo ',v6,v8,v9
             endif

 400       continue

c.....replace the a1 array with the a2 array
c.....thus eliminating the a1 array from memory,
c.....by working on only two longitude sections of data
c.....this program deos not require a large amount of
c.....memory space to run, although the output file may
c.....be rather large.
           do 410 la = 1,latot
             a1(la) = a2(la)
 410       continue
          endif
 100    continue

c.....normal level, no points added

      elseif ( level .eq. 2 ) then
        do 500 lo = 1, loend
          if ( lo .lt. lost ) then
c..... but do not process data if out of longitude range
            read(20,rec=lo)dummy
          elseif ( lo .eq. lost ) then
c.....read minimum longitude data and write vertices in
c.....latitude range
            read(20,rec=lo)(a1(la),la=1,latot)
            do 510 la = last, laend
              x = 60. * ( lo - 1)
              y = 60. * ( la - 1)
              write(21,*)'v ',x,y,a1(la)
 510        continue
          else
c.....read longitud data and write appropriate vertex data
            read(20,rec=lo)(a2(la),la=1,latot)
            do 520 la = last, laend
              x = 60. * ( lo -1)
              y = 60. * ( la -1)
              write(21,*)'v ',x,y,a2(la)
 520        continue
            do 530 la = last+1, laend
              e1 = a1( la -1)
              e2 = a2( la -1)
              e3 = a2(la)
              e4 = a1(la)
              v1 = (lo-lost-1) * (laend-last+1) + (la-last)
              v2 = v1 + (laend-last+1)
              v3 = v2 + 1
              v4 = v1 + 1
              if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
     &           e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
                write(21,*)'fo ',v1,v2,v3
                write(21,*)'fo ',v1,v3,v4
              else
                write(21,*)'fo ',v4,v1,v2
                write(21,*)'fo ',v4,v2,v3
              endif
 530        continue

          endif

 500    continue


c.....low resolution level
c.....some longitude data will be skipped depending on the
c.....skip factor entered at the start of the program. this
c.....results in less polygons but terrain resolution is lost.
c.....good for long range viewing of terrain when details
c.....would not be visible anyway

      elseif ( level .eq. 3 ) then

        do 600 lo = 1, loend
          if ( lo .lt. lost ) then
c.....do not process if longitude out of range
            read(20,rec=lo)dummy
          elseif ( lo .eq. lost ) then
c.....process data in range
c.....only process data in longitude and latitude which are
c.....integer skip factor multiples from the minimums.
            read(20,rec=lo)(a1(la),la=1,latot)
            do 610 la = last, laend-1, skip
              x = 60.0 * ( lo -1)
              y = 60.0 * ( la -1)
              write(21,*)'v ',x,y,a1(la)
 610        continue
            x = 60.0 * ( lo -1)
            y = 60.0 * ( laend -1)
            write(21,*)'v ',x,y,a1(laend)
          elseif(lo.gt.lost.and.lo.lt.loend)then
            if( mod((lo-lost),skip) .ne. 0 ) then
              read(20,rec=lo)dummy
            else
              read(20,rec=lo)(a2(la),la=1,latot)
              do 620 la = last, laend-1, skip
                if(mod(la-last,skip).eq.0)then
                  x = 60. * ( lo -1)
                  y = 60. * ( la -1)
                  write(21,*)'v ',x,y,a2(la)
                endif
 620          continue
              x = 60. * ( lo -1)
              y = 60. * ( laend -1)
              write(21,*)'v ',x,y,a2(laend)

c.....processing and creating polygons
              do 630 la = last+1, laend-1
                print*,'la =',la
                print*,'last =',last
                if( mod(la-last,skip).eq.0)then
                  e1 = a1(la-skip)
                  e2 = a2(la-skip)
                  e3 = a2(la)
                  e4 = a1(la)
                  print*,'e1 =',e1
                  print*,'e2 =',e2
                  print*,'e3 =',e3
                  print*,'e4 =',e4
                  if(mod(laend-last,skip).eq.0)then
                    inc = (laend-last)/skip + 1
                  else
                    inc = (laend-last)/skip + 2
                  endif
                  v1=((lo-lost)/skip-1)*inc + (la-last)/skip
                  v2 = v1 + inc
                  v3 = v2 + 1
                  v4 = v1 + 1
                  if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
     &               e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
                    write(21,*)'fo ',v1,v2,v3
                    write(21,*)'fo ',v1,v3,v4
                  else
                    write(21,*)'fo ',v4,v1,v2
                    write(21,*)'fo ',v4,v2,v3
                  endif
                endif
 630          continue
c.....compute polygons for la=laend value
c.....maximum and minimum values of latitude and longitude will
c.....be processed, even if the last geometric section will be
c.....a different size of the others, this allows the overall geometry
c.....to be the desired size with only small deviations appearing
c.....on, at most, two edges.
              e1=e4
              e2=e3
              e3=a2(laend)
              e4=a1(laend)
              print*,'e1 =',e1
              print*,'e2 =',e2
              print*,'e3 =',e3
              print*,'e4 =',e4
              v1=v4
              v2=v3
              v3=v2+1
              v4=v1+1
              print*,'v1 =',v1
              print*,'v2 =',v2
              print*,'v3 =',v3
              print*,'v4 =',v4
              if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
     &           e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
                write(21,*)'fo ',v1,v2,v3
                write(21,*)'fo ',v1,v3,v4
              else
                write(21,*)'fo ',v4,v1,v2
                write(21,*)'fo ',v4,v2,v3
              endif
c.....replace a1 array with a2 array
              do 640 la = 1, latot
                a1(la) = a2(la)
 640          continue
              losave = lo
              print*,'losave =',losave
            endif
          elseif ( lo .eq. loend ) then
c.....last value in range is always processed
            read(20,rec=lo)(a2(la),la=1,latot)
            do 660 la = last, laend-1, skip
              if(mod(la-last,skip).eq.0)then
                x = 60. * ( lo -1)
                y = 60. * ( la -1)
                write(21,*)'v ',x,y,a2(la)
              endif
 660        continue
            x = 60. * ( lo -1)
            y = 60. * ( laend -1)
            write(21,*)'v ',x,y,a2(laend)
              do 670 la = last+1, laend-1
                print*,'la =',la
                print*,'last =',last
                if( mod(la-last,skip).eq.0)then
                  e1 = a1(la-skip)
                  e2 = a2(la-skip)
                  e3 = a2(la)
                  e4 = a1(la)
                  print*,'e1 =',e1
                  print*,'e2 =',e2
                  print*,'e3 =',e3
                  print*,'e4 =',e4
                  if(mod(laend-last,skip).eq.0)then
                    inc = (laend-last)/skip + 1
                  else
                    inc = (laend-last)/skip + 2
                  endif
c...
                  v1=((losave-lost)/skip)*inc + (la-last)/skip
c...
                  v2 = v1 + inc
                  v3 = v2 + 1
                  v4 = v1 + 1
                  if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
     &               e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
                    write(21,*)'fo ',v1,v2,v3
                    write(21,*)'fo ',v1,v3,v4
                  else
                    write(21,*)'fo ',v4,v1,v2
                    write(21,*)'fo ',v4,v2,v3
                  endif
                endif
 670          continue
c.....compute polygons for la=laend value
              e1=e4
              e2=e3
              e3=a2(laend)
              e4=a1(laend)
              print*,'e1 =',e1
              print*,'e2 =',e2
              print*,'e3 =',e3
              print*,'e4 =',e4
              v1=v4
              v2=v3
              v3=v2+1
              v4=v1+1
              print*,'v1 =',v1
              print*,'v2 =',v2
              print*,'v3 =',v3
              print*,'v4 =',v4
              if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
     &           e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
                write(21,*)'fo ',v1,v2,v3
                write(21,*)'fo ',v1,v3,v4
              else
                write(21,*)'fo ',v4,v1,v2
                write(21,*)'fo ',v4,v2,v3
              endif
c.....replace a1 array with a2 array
              do 680 la = 1, latot
                a1(la) = a2(la)
 680          continue
        endif
 600    continue
      endif


 1000 end

c.....subroutine to compute the number of 3 second steps between
c.....two topographical longitude or latitude coordinates in the form
c..... ddd mm ss

      subroutine geodiff(max, min, istep)

      integer max(3), min(3)

c.....compute change in position

      if(max(3) .lt. min(3))then
        istep3 = 60 + max(3) - min(3)
        max(2) = max(2) - 1
      else
        istep3 = max(3) - min(3)
      endif

      if(max(2) .lt. min(2))then
        istep2 = 60 + max(2) - min(2)
        max(1) = max(1) - 1
      else
        istep2 = max(2) - min(2)
      endif

      istep1 = max(1) - min(1)

c......compute total number of steps

      istep = istep1*1200 + istep2*20 + istep3/3 + 1

      return

      end
