home *** CD-ROM | disk | FTP | other *** search
/ Avalon - 3D Objects & Resources / Avalon.iso / objects / obj / terrain / dmatoobj.f < prev    next >
Encoding:
Text File  |  1995-01-01  |  22.4 KB  |  677 lines
  1. c     *******************************************
  2. c     * Terrain Geometry Generation Program
  3. c     *
  4. c     * Kent F. Martin
  5. c     * 15 JAN 90
  6. c     * Animated Technologies
  7. c     * (213) 675 - 0770
  8. c     * 
  9. c     * This program reads Defense Terrain Mapping data and
  10. c     * converts is to a geometry file format which is compatible
  11. c     * with WaveFront Technology software (Model, Preview, Image)
  12. c     *
  13. c     * This version offers the choice of high, normal, or low
  14. c     * resolution output. Normal resolution output has polygons
  15. c     * with vertices located at the vertices in the database,
  16. c     * high resolution adds one vertex between the original
  17. c     * ones, resulting in 4 times more polygons, low resolution
  18. c     * skips some data values depending on user choice as to how 
  19. c     * much data to skip between saves, this results in considerably
  20. c     * less polygons
  21. c     *
  22. c     ******************************************
  23.  
  24. c.....Variable description list
  25.  
  26. c.....arrays:
  27. c       lomin : minimum longitude of database
  28. c       lomax : maximum longitude of database
  29. c       lamin : minimum latitude of database
  30. c       lamax : maximum latitude of database
  31. c
  32. c       lolo : minimum longitude of geometry region
  33. c       uplo : maximum longitude of geometry region
  34. c       lola : minimum latitude of geometry region
  35. c       upla : maximum latitude of geometry region
  36. c
  37. c       a1 : storage array #1, holds elevation data for processing
  38. c       a2 : storage array #2, holds elevation data for processing
  39. c            these are the only two storage arrays, making possible
  40. c            the processing of a large database without the need
  41. c            for large memory for program initialization.
  42. c
  43. c       losave : storage array for data when longitude values are being
  44. c                skipped in lo resolution, serves same purpose as a1 array
  45. c.....single variables
  46. c
  47. c     v1,v2,v3,...,v9 : temporary vertex numbers for polygon creation
  48. c
  49. c     xint : x coordinate midway between lo and (lo-1)
  50. c     yint : y coordinate midway between la and (la-1)
  51. c     int  : elevation at a point utilizing xint or yint, but not both
  52. c     int1 : elevation at a point when two intermediate points exist
  53. c     int2 : elevation at a point when two intermediate points exist
  54. c
  55. c
  56. c     lotot : counter of number of longitude values in database
  57. c     latot : number of latitude values in database
  58. c     lost : counter where processing of longitude data begins
  59. c     loend: counter where processing of longitude data ends
  60. c     last : counter where processing of latitude data begins
  61. c     laend: counter where processing of latitude data ends
  62.  
  63. c     ic   : integer seed value for intermediate point altitude
  64. c            generation
  65. c     z    : intermediate point altitude, result of using then
  66. c            Gaussian distribution random number generator
  67. c     ave  : average of elevations of points surrounding 
  68. c            intermediate point
  69. c     skip : skip factor for low resolution geometry
  70. c
  71. c     lo : longitude counter
  72. c     la : latitude counter
  73. c
  74. c.....file names
  75. c     in : file containing elevation data
  76. c     out: output file, this will be in a form which "model"
  77. c          can read, suggest user supply file.obj for name.
  78.  
  79.       character*80 in, out
  80.       real a1(5000), a2(5000)
  81.       real int, int1, int2
  82.       integer lolo(3), uplo(3), lola(3), upla(3)
  83.       integer lomax(3), lomin(3), lamax(3), lamin(3)
  84.       integer v1,v2,v3,v4,v5,v6,v7,v8,v9
  85.       integer skip
  86.       integer losave
  87.       integer choice
  88.  
  89.       print*,'Welcome to the Terrain Generator'
  90.       print*,' '
  91.       print*,'Enter data file name to be processed'
  92.       print*,' '
  93.       read(*,*)in
  94.       print*,'Enter output object file name'
  95.       print*,' '
  96.       read(*,*)out
  97.       print*,' '
  98.       print*,'Select desired level of terrain resolution:'
  99.       print*,'------'
  100.       print*,'  1) high'
  101.       print*,'  2) normal'
  102.       print*,'  3) low'
  103.       print*,' '
  104.       read(*,*)level
  105.  
  106. c.....if low level, ask for enlargement factor
  107.  
  108.       if ( level .eq. 3 ) then
  109.         print*,' '
  110.         print*,'Enter enlargement factor for geometry (integer)'
  111.         read(*,*)skip
  112.       endif
  113.  
  114. c.....read lower and upper limits for geometry boundary
  115.  
  116.       print*,'This program is intended to be used for elevation'
  117.       print*,'data recorded at 3 second intervals of longitude and'
  118.       print*,'latitude. If the value for seconds(ss) is not an'
  119.       print*,'integer multiple of 3, the input will be requested'
  120.       print*,'again.'
  121.  
  122. c.....enter database limits of longitude and latitude
  123.  6    print*,' '
  124.       print*,'Enter database lower longitude limit ( dd(d) mm ss )'
  125.       read(*,*)(lomin(i),i=1,3)
  126.       if(mod(lomin(3),3) .ne. 0 ) goto 6
  127.  
  128.  7    print*,' '
  129.       print*,'Enter database upper longitude limit ( dd(d) mm ss )'
  130.       read(*,*)(lomax(i),i=1,3)
  131.       if(mod(lomax(3),3) .ne. 0 ) goto 7
  132.  
  133.  8    print*,' '
  134.       print*,'Enter database lower latitude limit ( dd mm ss )'
  135.       read(*,*)(lamin(i),i=1,3)
  136.       if(mod(lamin(3),3) .ne. 0 ) goto 8
  137.  
  138.  9    print*,' '
  139.       print*,'Enter database upper latitude limit ( dd mm ss )'
  140.       read(*,*)(lamax(i),i=1,3)
  141.       if(mod(lamax(3),3) .ne. 0 ) goto 9
  142.  
  143. c.....specify area of database to be processed
  144.  
  145.       print*,' '
  146.  10   print*,'Enter lower longitude limit ( dd(d) mm ss )'
  147.       read(*,*)(lolo(i),i=1,3)
  148.       if( mod(lolo(3),3) .ne. 0 )goto 10
  149.  
  150.       print*,' '
  151.  11   print*,'Enter upper longitude limit ( dd(d) mm ss )'
  152.       read(*,*)(uplo(i),i=1,3)
  153.       if( mod(uplo(3),3) .ne. 0 )goto 11
  154.  
  155.       print*,' '
  156.  12   print*,'Enter lower latitude limit ( dd mm ss )'
  157.       read(*,*)(lola(i),i=1,3)
  158.       if( mod(lola(3),3) .ne. 0)goto 12
  159.  
  160.       print*,' '
  161.  13   print*,'Enter upper latitude limit ( dd mm ss )'
  162.       read(*,*)(upla(i),i=1,3)
  163.       if( mod(upla(3),3) .ne. 0)goto13
  164.  
  165. c....., form='unformated'open input and output files
  166.  
  167.       open(unit=20, file=in, status='old',form='unformatted',
  168.      & access='direct',recl=9604)
  169.       open(unit=21, file=out, status='unknown')
  170.  
  171. c.....establish do loop control counters by using geodiff subroutine
  172. c.....minimum and maximum longitude and latitude should be read in from
  173. c.....data file header, if not it must be requested interactively
  174. c.....they will be hard coded for this application
  175.  
  176.       call geodiff(lolo, lomin, lost)
  177.       call geodiff(uplo, lomin, loend)
  178.       call geodiff(lola, lamin, last)
  179.       call geodiff(upla, lamin, laend)
  180.       call geodiff(lomax, lomin, lotot)
  181.       call geodiff(lamax, lamin, latot)
  182. c.....print control counters for debugging
  183. c     print*,'lotot =',lotot
  184. c     print*,'latot =',latot
  185. c     print*,'lost =',lost
  186. c     print*,'loend =',loend
  187. c     print*,'last =',last
  188. c     print*,'laend =',laend
  189.  
  190. c.....compute number of vertices and polygons
  191. c.....notify user of these
  192. c.....give option to reenter limits, continue, or quit.
  193.  
  194.       if(level.eq.1)then
  195.         numv=((loend-lost+1)*2-1)*((laend-last+1)*2-1)
  196.         nump=(loend-lost)*(laend-last)*8
  197.       elseif(level.eq.2)then
  198.         numv=(loend-lost+1)*(laend-last+1)
  199.         nump=(loend-lost)*(laend-last)*2
  200.       elseif(level.eq.3)then
  201.         if(mod(loend-lost,skip) .eq. 0 )then
  202.           loadd = 1
  203.         else
  204.           loadd = 2
  205.         endif
  206.         if(mod(laend-last,skip) .eq. 0 )then
  207.           laadd = 1
  208.         else
  209.           laadd = 2
  210.         endif
  211.         numv=((loend-lost)/skip+loadd)*((laend-last)/skip+laadd)
  212.         nump=((loend-lost)/skip+loadd-1)*
  213.      &       ((laend-last) /skip+laadd-1) * 2
  214.       endif
  215.  
  216.       print*,' '
  217.       print*,'Current geography limits and resolution type'
  218.       print*,'        will yield:'
  219.       print*,' ',numv,' vertices and'
  220.       print*,' ',nump,' polygons.'
  221.       print*,' '
  222.  80   print*,'Do you wish to:'
  223.       print*,'  1) continue'
  224.       print*,'  2) enter new geography limits'
  225.       print*,'  3) quit?'
  226.       read(*,*)choice
  227.  
  228.       if(choice .eq. 1)then
  229.         goto 90
  230.       elseif(choice .eq. 2)then
  231.         goto 6
  232.       elseif(choice .eq. 3)then
  233.         goto 1000
  234.       else
  235.         goto 80
  236.       endif
  237.  
  238. c.....branch to correct area depending on resolution level
  239.  
  240. c.....high resolution area
  241.  
  242.  90   continue
  243.       if ( level .eq. 1 ) then
  244.  
  245. c.....do only enough lines to cover longitude range
  246.  
  247.         do 100 lo = 1, loend
  248.  
  249.           if ( lo .lt. lost ) then
  250. c..... a dummy line if longitude out of range
  251.             read(20,rec=lo)dummy
  252. c.....            format(a1)
  253.           elseif ( lo .eq. lost ) then
  254. c.....read first line of data and write vertex information
  255.             read(20,rec=lo)(a1(la),la=1,latot)
  256.             do 110 la = last, laend
  257.               if ( la .eq. last ) then
  258. c.....write only one point if latitude is minimum value
  259.                 x = 60. * ( lo - 1)
  260.                 y = 60. * ( la - 1)
  261.                 write(21,*)'v ',x,y,a1(la)
  262.               else
  263. c.....write two point if latitude is greater than minimum
  264.                 x = 60. * ( lo -1)
  265.                 y = 60. * ( la -1)
  266.                 yint = y - 30.
  267. c.....intermediate grid point #1
  268.                 ic = lo * ( la - 1 )
  269.                 z = grn(ic)
  270.                 ave = (a1(la) + a1(la-1)) / 2.
  271.                 deltaz = abs( a1(la) - a1(la-1) )
  272.                 int = ave + deltaz / 10. * z
  273.                 write(21,*)'v ',x,yint,int
  274.                 write(21,*)'v ',x,y,a1(la)
  275.               endif
  276.  110        continue
  277.           else
  278. c.....if longitude is not starting value, ad a line of data
  279. c.....and store in second array
  280.            read(20,rec=lo)(a2(la),la=1,latot)
  281.            do 200 la = last, laend
  282.              if ( la .eq. last ) then
  283. c.....write two vertices for minimum latitude
  284.                xint = 60.0 * ( lo - 1) - 30.
  285.                y = 60. * (la -1)
  286. c.....intermediate grid point #2
  287.                ic = (lo - 1) * ( la - 1)
  288.                z = grn(ic)
  289.                ave = ( a1(la) + a2(la) ) / 2.
  290.                deltaz = abs ( a1(la) - a2(la) )
  291.                int = ave + deltaz / 10. * z
  292.                write(21,*)'v ',xint,y,int
  293.              else
  294. c.....determine and write two more intermediate points
  295. c.....intermediate grid point #3
  296.                ic = (lo-1)*(la-1)+(lo*la)
  297.                z = grn(ic)
  298.                ave =(a1(la)+a2(la)+a1(la-1)+a2(la-1))/4.
  299.                zmax = max(a1(la),a2(la),a1(la-1),a2(la-1))
  300.                zmin = min(a1(la),a2(la),a1(la-1),a2(la-1))
  301.                deltaz = abs(zmax-zmin)
  302.                int1 = ave + deltaz / 10. * z
  303. c.....intermediate grid point #4
  304.                ic = (lo-1)*la
  305.                z = grn(ic)
  306.                ave = (a1(la) + a2(la)) / 2.
  307.                deltaz = abs( a1(la) - a2(la) )
  308.                int2 = ave + deltaz / 10. * z
  309.                xint = 60. * ( lo -1) - 30.
  310.                yint = 60. * ( la -1) - 30.
  311.                y = 60. * (la -1)
  312.                write(21,*)'v ',xint,yint,int1
  313.                write(21,*)'v ',xint,y,int2
  314.              endif
  315.  200       continue
  316.  
  317. c.....without increasin longitude counter, go through latitude
  318. c.....values generating the right hand side of the geometry 
  319. c.....generated to this point.
  320.            do 300 la = last, laend
  321.              if ( la .eq. last ) then
  322.                x = 60. * ( lo -1)
  323.                y = 60. * ( la -1)
  324.                write(21,*)'v ',x,y,a2(la)
  325.              else
  326.                x = 60. * ( lo -1)
  327.                y = 60. * ( la -1)
  328.                yint = y - 30.
  329. c.....intermediate point #5
  330.                ic = lo * la
  331.                z = grn(ic)
  332.                ave = ( a2(la) + a2(la-1) ) / 2.
  333.                deltaz = abs ( a2(la) - a2(la-1) )
  334.                int = ave + deltaz / 10. * z
  335.                write(21,*)'v ',x,yint,int
  336.                write(21,*)'v ',x,y,a2(la)
  337.              endif
  338.  300       continue
  339.  
  340. c.....begin the creation of polygons
  341. c.....start with a longitude value one step larger than minimum
  342.  
  343. c.....polygons are built by analyzing four points in the
  344. c.....original database which make a rectangle. this region
  345. c.....is converted into 8 triangles by use of the intermediate
  346. c.....points generated earlier. the 'direction' of the triangles
  347. c.....is determined by finding the highest corner and working
  348. c.....down from that point.
  349.            do 400 la = last+1, laend
  350.              e1 = a1(la-1)
  351.              e2 = a2(la-1)
  352.              e3 = a2(la)
  353.              e4 = a1(la)
  354.              print*,'lo =',lo
  355.              print*,'la =',la
  356.              print*,'lost =',lost
  357.              print*,'last =',last
  358.              print*,'laend =',laend
  359.              v1=(lo-lost-1)*(((laend-last+1)*2-1)*2)+(la-last)*2 - 1
  360.              v2 = v1 + 1
  361.              v3 = v1 + 2
  362.              v4 = v1 + (((laend-last+1)*2)-1)
  363.              v5 = v4 + 1
  364.              v6 = v4 + 2
  365.              v7 = v4 + (((laend-last+1)*2)-1)
  366.              v8 = v7 + 1
  367.              v9 = v7 + 2
  368.              if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
  369.      &          e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
  370.                write(21,*)'fo ',v1,v4,v5
  371.                write(21,*)'fo ',v1,v5,v2
  372.                write(21,*)'fo ',v2,v5,v6
  373.                write(21,*)'fo ',v2,v6,v3
  374.                write(21,*)'fo ',v4,v7,v8
  375.                write(21,*)'fo ',v4,v8,v5
  376.                write(21,*)'fo ',v5,v8,v9
  377.                write(21,*)'fo ',v5,v9,v6
  378.              else
  379.                write(21,*)'fo ',v2,v1,v4
  380.                write(21,*)'fo ',v2,v4,v5
  381.                write(21,*)'fo ',v3,v2,v5
  382.                write(21,*)'fo ',v3,v5,v6
  383.                write(21,*)'fo ',v5,v4,v7
  384.                write(21,*)'fo ',v5,v7,v8
  385.                write(21,*)'fo ',v6,v5,v8
  386.                write(21,*)'fo ',v6,v8,v9
  387.              endif
  388.  
  389.  400       continue
  390.  
  391. c.....replace the a1 array with the a2 array
  392. c.....thus eliminating the a1 array from memory,
  393. c.....by working on only two longitude sections of data
  394. c.....this program deos not require a large amount of
  395. c.....memory space to run, although the output file may
  396. c.....be rather large.
  397.            do 410 la = 1,latot
  398.              a1(la) = a2(la)
  399.  410       continue
  400.           endif
  401.  100    continue
  402.  
  403. c.....normal level, no points added
  404.  
  405.       elseif ( level .eq. 2 ) then
  406.         do 500 lo = 1, loend
  407.           if ( lo .lt. lost ) then
  408. c..... but do not process data if out of longitude range
  409.             read(20,rec=lo)dummy
  410.           elseif ( lo .eq. lost ) then
  411. c.....read minimum longitude data and write vertices in
  412. c.....latitude range
  413.             read(20,rec=lo)(a1(la),la=1,latot)
  414.             do 510 la = last, laend
  415.               x = 60. * ( lo - 1)
  416.               y = 60. * ( la - 1)
  417.               write(21,*)'v ',x,y,a1(la)
  418.  510        continue
  419.           else
  420. c.....read longitud data and write appropriate vertex data
  421.             read(20,rec=lo)(a2(la),la=1,latot)
  422.             do 520 la = last, laend
  423.               x = 60. * ( lo -1)
  424.               y = 60. * ( la -1)
  425.               write(21,*)'v ',x,y,a2(la)
  426.  520        continue
  427.             do 530 la = last+1, laend
  428.               e1 = a1( la -1)
  429.               e2 = a2( la -1)
  430.               e3 = a2(la)
  431.               e4 = a1(la)
  432.               v1 = (lo-lost-1) * (laend-last+1) + (la-last)
  433.               v2 = v1 + (laend-last+1)
  434.               v3 = v2 + 1
  435.               v4 = v1 + 1
  436.               if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
  437.      &           e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
  438.                 write(21,*)'fo ',v1,v2,v3
  439.                 write(21,*)'fo ',v1,v3,v4
  440.               else
  441.                 write(21,*)'fo ',v4,v1,v2
  442.                 write(21,*)'fo ',v4,v2,v3
  443.               endif
  444.  530        continue
  445.  
  446.           endif
  447.  
  448.  500    continue
  449.  
  450.  
  451. c.....low resolution level
  452. c.....some longitude data will be skipped depending on the
  453. c.....skip factor entered at the start of the program. this
  454. c.....results in less polygons but terrain resolution is lost.
  455. c.....good for long range viewing of terrain when details
  456. c.....would not be visible anyway
  457.  
  458.       elseif ( level .eq. 3 ) then
  459.  
  460.         do 600 lo = 1, loend
  461.           if ( lo .lt. lost ) then
  462. c.....do not process if longitude out of range
  463.             read(20,rec=lo)dummy
  464.           elseif ( lo .eq. lost ) then
  465. c.....process data in range
  466. c.....only process data in longitude and latitude which are
  467. c.....integer skip factor multiples from the minimums.
  468.             read(20,rec=lo)(a1(la),la=1,latot)
  469.             do 610 la = last, laend-1, skip
  470.               x = 60.0 * ( lo -1)
  471.               y = 60.0 * ( la -1)
  472.               write(21,*)'v ',x,y,a1(la)
  473.  610        continue
  474.             x = 60.0 * ( lo -1)
  475.             y = 60.0 * ( laend -1)
  476.             write(21,*)'v ',x,y,a1(laend)
  477.           elseif(lo.gt.lost.and.lo.lt.loend)then
  478.             if( mod((lo-lost),skip) .ne. 0 ) then
  479.               read(20,rec=lo)dummy
  480.             else
  481.               read(20,rec=lo)(a2(la),la=1,latot)
  482.               do 620 la = last, laend-1, skip
  483.                 if(mod(la-last,skip).eq.0)then
  484.                   x = 60. * ( lo -1)
  485.                   y = 60. * ( la -1)
  486.                   write(21,*)'v ',x,y,a2(la)
  487.                 endif
  488.  620          continue
  489.               x = 60. * ( lo -1)
  490.               y = 60. * ( laend -1)
  491.               write(21,*)'v ',x,y,a2(laend)
  492.  
  493. c.....processing and creating polygons
  494.               do 630 la = last+1, laend-1
  495.                 print*,'la =',la
  496.                 print*,'last =',last
  497.                 if( mod(la-last,skip).eq.0)then
  498.                   e1 = a1(la-skip)
  499.                   e2 = a2(la-skip)
  500.                   e3 = a2(la)
  501.                   e4 = a1(la)
  502.                   print*,'e1 =',e1
  503.                   print*,'e2 =',e2
  504.                   print*,'e3 =',e3
  505.                   print*,'e4 =',e4
  506.                   if(mod(laend-last,skip).eq.0)then
  507.                     inc = (laend-last)/skip + 1
  508.                   else
  509.                     inc = (laend-last)/skip + 2
  510.                   endif
  511.                   v1=((lo-lost)/skip-1)*inc + (la-last)/skip
  512.                   v2 = v1 + inc
  513.                   v3 = v2 + 1
  514.                   v4 = v1 + 1
  515.                   if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
  516.      &               e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
  517.                     write(21,*)'fo ',v1,v2,v3
  518.                     write(21,*)'fo ',v1,v3,v4
  519.                   else
  520.                     write(21,*)'fo ',v4,v1,v2
  521.                     write(21,*)'fo ',v4,v2,v3
  522.                   endif
  523.                 endif
  524.  630          continue
  525. c.....compute polygons for la=laend value
  526. c.....maximum and minimum values of latitude and longitude will
  527. c.....be processed, even if the last geometric section will be
  528. c.....a different size of the others, this allows the overall geometry
  529. c.....to be the desired size with only small deviations appearing
  530. c.....on, at most, two edges.
  531.               e1=e4
  532.               e2=e3
  533.               e3=a2(laend)
  534.               e4=a1(laend)
  535.               print*,'e1 =',e1
  536.               print*,'e2 =',e2
  537.               print*,'e3 =',e3
  538.               print*,'e4 =',e4
  539.               v1=v4
  540.               v2=v3
  541.               v3=v2+1
  542.               v4=v1+1
  543.               print*,'v1 =',v1
  544.               print*,'v2 =',v2
  545.               print*,'v3 =',v3
  546.               print*,'v4 =',v4
  547.               if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
  548.      &           e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
  549.                 write(21,*)'fo ',v1,v2,v3
  550.                 write(21,*)'fo ',v1,v3,v4
  551.               else
  552.                 write(21,*)'fo ',v4,v1,v2
  553.                 write(21,*)'fo ',v4,v2,v3
  554.               endif
  555. c.....replace a1 array with a2 array
  556.               do 640 la = 1, latot
  557.                 a1(la) = a2(la)
  558.  640          continue
  559.               losave = lo
  560.               print*,'losave =',losave
  561.             endif
  562.           elseif ( lo .eq. loend ) then
  563. c.....last value in range is always processed
  564.             read(20,rec=lo)(a2(la),la=1,latot)
  565.             do 660 la = last, laend-1, skip
  566.               if(mod(la-last,skip).eq.0)then
  567.                 x = 60. * ( lo -1)
  568.                 y = 60. * ( la -1)
  569.                 write(21,*)'v ',x,y,a2(la)
  570.               endif
  571.  660        continue
  572.             x = 60. * ( lo -1)
  573.             y = 60. * ( laend -1)
  574.             write(21,*)'v ',x,y,a2(laend)
  575.               do 670 la = last+1, laend-1
  576.                 print*,'la =',la
  577.                 print*,'last =',last
  578.                 if( mod(la-last,skip).eq.0)then
  579.                   e1 = a1(la-skip)
  580.                   e2 = a2(la-skip)
  581.                   e3 = a2(la)
  582.                   e4 = a1(la)
  583.                   print*,'e1 =',e1
  584.                   print*,'e2 =',e2
  585.                   print*,'e3 =',e3
  586.                   print*,'e4 =',e4
  587.                   if(mod(laend-last,skip).eq.0)then
  588.                     inc = (laend-last)/skip + 1
  589.                   else
  590.                     inc = (laend-last)/skip + 2
  591.                   endif
  592. c...
  593.                   v1=((losave-lost)/skip)*inc + (la-last)/skip
  594. c...
  595.                   v2 = v1 + inc
  596.                   v3 = v2 + 1
  597.                   v4 = v1 + 1
  598.                   if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
  599.      &               e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
  600.                     write(21,*)'fo ',v1,v2,v3
  601.                     write(21,*)'fo ',v1,v3,v4
  602.                   else
  603.                     write(21,*)'fo ',v4,v1,v2
  604.                     write(21,*)'fo ',v4,v2,v3
  605.                   endif
  606.                 endif
  607.  670          continue
  608. c.....compute polygons for la=laend value
  609.               e1=e4
  610.               e2=e3
  611.               e3=a2(laend)
  612.               e4=a1(laend)
  613.               print*,'e1 =',e1
  614.               print*,'e2 =',e2
  615.               print*,'e3 =',e3
  616.               print*,'e4 =',e4
  617.               v1=v4
  618.               v2=v3
  619.               v3=v2+1
  620.               v4=v1+1
  621.               print*,'v1 =',v1
  622.               print*,'v2 =',v2
  623.               print*,'v3 =',v3
  624.               print*,'v4 =',v4
  625.               if(e1.ge.e2.and.e1.ge.e3.and.e1.ge.e4.or.
  626.      &           e3.ge.e1.and.e3.ge.e2.and.e3.ge.e4)then
  627.                 write(21,*)'fo ',v1,v2,v3
  628.                 write(21,*)'fo ',v1,v3,v4
  629.               else
  630.                 write(21,*)'fo ',v4,v1,v2
  631.                 write(21,*)'fo ',v4,v2,v3
  632.               endif
  633. c.....replace a1 array with a2 array
  634.               do 680 la = 1, latot
  635.                 a1(la) = a2(la)
  636.  680          continue
  637.         endif
  638.  600    continue
  639.       endif
  640.  
  641.  
  642.  1000 end
  643.  
  644. c.....subroutine to compute the number of 3 second steps between
  645. c.....two topographical longitude or latitude coordinates in the form
  646. c..... ddd mm ss
  647.  
  648.       subroutine geodiff(max, min, istep)
  649.  
  650.       integer max(3), min(3)
  651.  
  652. c.....compute change in position
  653.  
  654.       if(max(3) .lt. min(3))then
  655.         istep3 = 60 + max(3) - min(3)
  656.         max(2) = max(2) - 1
  657.       else
  658.         istep3 = max(3) - min(3)
  659.       endif
  660.  
  661.       if(max(2) .lt. min(2))then
  662.         istep2 = 60 + max(2) - min(2)
  663.         max(1) = max(1) - 1
  664.       else
  665.         istep2 = max(2) - min(2)
  666.       endif
  667.  
  668.       istep1 = max(1) - min(1)
  669.  
  670. c......compute total number of steps
  671.  
  672.       istep = istep1*1200 + istep2*20 + istep3/3 + 1
  673.  
  674.       return
  675.  
  676.       end
  677.