1cb::inverse
2c
3      program inverse
4c
5c********1*********2*********3*********4*********5*********6*********7**
6c
7c name:      inverse
8c version:   201211.29
9c author:    stephen j. frakes
10c last mod:  Charles Karney
11c purpose:   to compute a geodetic inverse
12c            and then display output information
13c
14c input parameters:
15c -----------------
16c
17c output parameters:
18c ------------------
19c
20c local variables and constants:
21c ------------------------------
22c answer           user prompt response
23c b                semiminor axis polar (in meters)
24c baz              azimuth back (in radians)
25c buff             input char buffer array
26c dd,dm,ds         temporary values for degrees, minutes, seconds
27c dlon             temporary value for difference in longitude (radians)
28c dmt,d_mt         char constants for meter units
29c edist            ellipsoid distance (in meters)
30c elips            ellipsoid choice
31c esq              eccentricity squared for reference ellipsoid
32c faz              azimuth forward (in radians)
33c filout           output file name
34c finv             reciprocal flattening
35c hem              hemisphere flag for lat & lon entry
36c ierror           error condition flag with d,m,s conversion
37c lgh              length of buff() array
38c option           user prompt response
39c r1,r2            temporary variables
40c ss               temporary variable
41c tol              tolerance for conversion of seconds
42c
43c name1            name of station one
44c ld1,lm1,sl1      latitude  sta one - degrees,minutes,seconds
45c ald1,alm1,sl1    latitude  sta one - degrees,minutes,seconds
46c lat1sn           latitude  sta one - sign (+/- 1)
47c d_ns1            latitude  sta one - char ('N','S')
48c lod1,lom1,slo1   longitude sta one - degrees,minutes,seconds
49c alod1,alom1,slo1 longitude sta one - degrees,minutes,seconds
50c lon1sn           longitude sta one - sign (+/- 1)
51c d_ew1            longitude sta one - char ('E','W')
52c iaz1,maz1,saz1   forward azimuth   - degrees,minutes,seconds
53c isign1           forward azimuth   - flag  (+/- 1)
54c glat1,glon1      station one       - (lat & lon in radians )
55c p1,e1            standpoint one    - (lat & lon in radians )
56c
57c name2            name of station two
58c ld2,lm2,sl2      latitude  sta two - degrees,minutes,seconds
59c ald2,alm2,sl2    latitude  sta two - degrees,minutes,seconds
60c lat2sn           latitude  sta two - sign (+/- 1)
61c d_ns2            latitude  sta two - char ('N','S')
62c lod2,lom2,slo2   longitude sta two - degrees,minutes,seconds
63c alod2,alom2,slo2 longitude sta one - degrees,minutes,seconds
64c lon2sn           longitude sta two - sign (+/- 1)
65c d_ew2            longitude sta two - char ('E','W')
66c iaz2,maz2,saz2   back azimuth      - degrees,minutes,seconds
67c isign2           back azimuth      - flag  (+/- 1)
68c glat2,glon2      station two       - (lat & lon in radians )
69c p2,e2            forepoint two     - (lat & lon in radians )
70c
71c global variables and constants:
72c -------------------------------
73c a                semimajor axis equatorial (in meters)
74c f                flattening
75c pi               constant 3.14159....
76c rad              constant 180.0/pi
77c
78c    this module called by:  n/a
79c
80c    this module calls:      elipss, getdeg, inver1, todmsp
81c    gethem, trim,   bufdms, gvalr8, gvali4, fixdms, gpnhri ***********
82c    gethem, trim,   bufdms, gvalr8, gvali4, fixdms, invers <----------
83c    datan,  write,  read,   dabs,   open,   stop
84c
85c    include files used:     n/a
86c
87c    common blocks used:     const, elipsoid
88c
89c    references:             see comments within subroutines
90c
91c    comments:
92c
93c********1*********2*********3*********4*********5*********6*********7**
94c::modification history
95c::1990mm.dd, sjf, creation of program
96c::199412.15, bmt, creation of program on viper
97c::200203.08, crs, modified by c.schwarz to correct spelling of Clarke
98c::200207.15, rws, modified i/o & standardized program documentation
99c::                added subs trim, bufdms, gethem, gvali4, gvalr8
100c::200207.23, rws, replaced sub inver1 with gpnarc, gpnloa, gpnhri
101c::200208.15, rws, fixed an error in bufdms
102c::              - renamed ellips to elipss "common error" with dirct2
103c::              - added FAZ & BAZ to printed output
104c::200208.19, rws, added more error flags for web interface code
105c::              - added logical nowebb
106c::200208.xx, sjf, program version number 2.0
107c::201105.xx, dgm, program version number 3.0
108c::              - replaced sub gpnarc, gpnloa, gpnhri with invers
109c::              - tested for valid antipodal solutions (+/- 0.1 mm)
110c::              - tested for polar solutions (+/- 0.1 mm)
111c::              - needs improvement for long-line/antipodal boundary
112c::201211.29, cffk, program version numer 3.1
113c::              - drop in replacement routines from
114c::                "Algorithms for Geodesics"
115c********1*********2*********3*********4*********5*********6*********7**
116ce::inverse
117c
118      implicit double precision (a-h, o-z)
119      implicit integer (i-n)
120c
121      logical  nowebb
122c
123      character*1  answer,option,dmt,buff(50),hem
124      character*6  d_ns1, d_ew1, d_ns2, d_ew2, d_mt
125      character*30 filout,name1,name2,elips
126c
127      integer*4    ierror
128      integer*4    lgh
129c
130      common/const/pi,rad
131      common/elipsoid/a,f
132c
133c     ms_unix      0 = web version
134c                  1 = ms_dos or unix version
135c
136      parameter   ( ms_unix = 0 )
137c
138      pi   = 4.d0*datan(1.d0)
139      rad  = 180.d0/pi
140      dmt  = 'm'
141      d_mt = 'Meters'
142c
143      if( ms_unix.eq.1 )then
144        nowebb = .true.
145      else
146        nowebb = .false.
147      endif
148c
149    1 do 2 i=1,25
150        write(*,*) '  '
151    2 continue
152c
153    5 write(*,*) '  Program Inverse  -  Version 3.1 '
154      write(*,*) '  '
155      write(*,*) '  Ellipsoid options : '
156      write(*,*) '  '
157      write(*,*) '  1) GRS80 / WGS84  (NAD83) '
158      write(*,*) '  2) Clarke 1866    (NAD27) '
159      write(*,*) '  3) Any other ellipsoid '
160      write(*,*) '  '
161      write(*,*) '  Enter choice : '
162      read(*,10) option
163   10 format(a1)
164c
165      if(option.eq.'1') then
166        a=6378137.d0
167        f=1.d0/298.257222100882711243162836600094d0
168        elips='GRS80 / WGS84  (NAD83)'
169      elseif(option.eq.'2') then
170        a=6378206.4d0
171        f=1.d0/294.9786982138d0
172        elips='Clarke 1866    (NAD27)'
173      elseif(option.eq.'3') then
174        call elipss (elips)
175      else
176        write(*,*) '  Enter 1, 2, or 3 !   Try again --'
177        goto 5
178      endif
179c
180      esq = f*(2.0d0-f)
181c
182      write(*,*) '  '
183      write(*,*) '  '
184      write(*,*) '  '
185      write(*,*) '  '
186c
187   15 write(*,*) '  Enter First Station '
188      write(*,16)
189   16 format(18x,'(Separate D,M,S by blanks or commas)')
190      write(*,*) 'hDD MM SS.sssss  Latitude :        (h default = N )'
191   11 format(50a1)
192c
193   22 hem='N'
194      read(*,11) buff
195      call trim (buff,lgh,hem)
196      call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
197c
198      if( ierror.eq.0 )then
199        irlat1 = 0
200      else
201        irlat1 = 1
202        write(*,*) ' Invalid Latitude ... Please re-enter it '
203        write(*,*) '  '
204        if( nowebb )then
205          goto 22
206        endif
207      endif
208c
209      ald1 = dd
210      alm1 = dm
211      sl1  = ds
212c
213      if( hem.eq.'N' ) then
214        lat1sn = +1
215      else
216        lat1sn = -1
217      endif
218c
219      write(*,*) 'hDDD MM SS.sssss Longitude :       (h default = W )'
220c
221   23 hem='W'
222      read(*,11) buff
223      call trim (buff,lgh,hem)
224      call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
225c
226      if( ierror.eq.0 )then
227        irlon1 = 0
228      else
229        irlon1 = 1
230        write(*,*) ' Invalid Longitude ... Please re-enter it '
231        write(*,*) '  '
232        if( nowebb )then
233          goto 23
234        endif
235      endif
236c
237      alod1 = dd
238      alom1 = dm
239      slo1  = ds
240c
241      if( hem.eq.'E' ) then
242        lon1sn = +1
243      else
244        lon1sn = -1
245      endif
246c
247      call getdeg(ald1, alm1, sl1, lat1sn, glat1)
248      call getdeg(alod1,alom1,slo1,lon1sn, glon1)
249c
250      write(*,*) '  '
251      write(*,*) '  '
252c
253   20 write(*,*) '  Enter Second Station '
254      write(*,16)
255      write(*,*) 'hDD MM SS.sssss  Latitude :        (h default = N )'
256c
257   24 hem='N'
258      read(*,11) buff
259      call trim (buff,lgh,hem)
260      call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
261c
262      if( ierror.eq.0 )then
263        irlat2 = 0
264      else
265        irlat2 = 1
266        write(*,*) ' Invalid Latitude ... Please re-enter it '
267        write(*,*) '  '
268        if( nowebb )then
269          goto 24
270        endif
271      endif
272c
273      ald2 = dd
274      alm2 = dm
275      sl2  = ds
276c
277      if( hem.eq.'N' ) then
278        lat2sn = +1
279      else
280        lat2sn = -1
281      endif
282c
283      write(*,*) 'hDDD MM SS.sssss Longitude :       (h default = W )'
284c
285   25 hem='W'
286      read(*,11) buff
287      call trim (buff,lgh,hem)
288      call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
289c
290      if( ierror.eq.0 )then
291        irlon2 = 0
292      else
293        irlon2 = 1
294        write(*,*) ' Invalid Longitude ... Please re-enter it '
295        write(*,*) '  '
296        if( nowebb )then
297          goto 25
298        endif
299      endif
300c
301      alod2 = dd
302      alom2 = dm
303      slo2  = ds
304c
305      if( hem.eq.'E' )then
306        lon2sn = +1
307      else
308        lon2sn = -1
309      endif
310c
311      call getdeg(ald2, alm2, sl2, lat2sn, glat2)
312      call getdeg(alod2,alom2,slo2,lon2sn, glon2)
313c
314      p1 = glat1
315      e1 = glon1
316      p2 = glat2
317      e2 = glon2
318c
319      if( e1.lt.0.0d0 )then
320        e1 = e1+2.0d0*pi
321      endif
322      if( e2.lt.0.0d0 )then
323        e2 = e2+2.0d0*pi
324      endif
325c
326c     compute the geodetic inverse
327c
328      call invers(a, f, p1, e1, p2, e2,
329     +    edist, faz, baz, 0, dummy, dummy, dummy, dummy, dummy)
330      if (baz .ge. 0) then
331        baz = baz - 180
332      else
333        baz = baz + 180
334      end if
335c
336c     set the tolerance (in seconds) for the azimuth conversion
337c
338      tol = 0.00005d0
339c
340      call todmsp(faz,iaz1,maz1,saz1,isign1)
341      if(isign1.lt.0) then
342        iaz1=359-iaz1
343        maz1=59-maz1
344        saz1=60.d0-saz1
345      endif
346      call fixdms ( iaz1, maz1, saz1, tol )
347c
348      call todmsp(baz,iaz2,maz2,saz2,isign2)
349      if(isign2.lt.0) then
350        iaz2=359-iaz2
351        maz2=59-maz2
352        saz2=60.d0-saz2
353      endif
354      call fixdms ( iaz2, maz2, saz2, tol )
355c
356      call todmsp(glat1, ld1,  lm1,  sl1,  lat1sn)
357      call todmsp(glon1, lod1, lom1, slo1, lon1sn)
358      call todmsp(glat2, ld2,  lm2,  sl2,  lat2sn)
359      call todmsp(glon2, lod2, lom2, slo2, lon2sn)
360c
361      call hem_ns ( lat1sn, d_ns1 )
362      call hem_ew ( lon1sn, d_ew1 )
363      call hem_ns ( lat2sn, d_ns2 )
364      call hem_ew ( lon2sn, d_ew2 )
365c
366      write(*,*) '  '
367      write(*,*) '  '
368      write(*,*) '  '
369      write(*,*) '  '
370      write(*,*) '  '
371      write(*,49) elips
372   49 format('  Ellipsoid : ',a30)
373      finv=1.d0/f
374      b=a*(1.d0-f)
375      write(*,50) a,b,finv
376   50 format('  Equatorial axis,    a   = ',f15.4,/,
377     *       '  Polar axis,         b   = ',f15.4,/,
378     *       '  Inverse flattening, 1/f = ',f16.11)
379c
380   18 format('    LAT = ',i3,1x,i2,1x,f8.5,1x,a6)
381   19 format('    LON = ',i3,1x,i2,1x,f8.5,1x,a6)
382c
383      write(*,*) '  '
384      write(*,*) '  First  Station : '
385      write(*,*) '  ---------------- '
386      write(*,18) ld1, lm1, sl1, d_ns1
387      write(*,19) lod1,lom1,slo1,d_ew1
388c
389      write(*,*) '  '
390      write(*,*) '  Second Station : '
391      write(*,*) '  ---------------- '
392      write(*,18) ld2, lm2, sl2, d_ns2
393      write(*,19) lod2,lom2,slo2,d_ew2
394c
395   32 format('  Ellipsoidal distance     S = ',f14.4,1x,a1)
396   34 format('  Forward azimuth        FAZ = ',i3,1x,i2,1x,f7.4,
397     1       ' From North')
398   35 format('  Back azimuth           BAZ = ',i3,1x,i2,1x,f7.4,
399     1       ' From North')
400c
401      write(*,*) '  '
402      write(*,34) iaz1,maz1,saz1
403      write(*,35) iaz2,maz2,saz2
404      write(*,32) edist,dmt
405      write(*,*) '  '
406      write(*,*) '  Do you want to save this output into a file (y/n)?'
407      read(*,10) answer
408c
409      if( answer.eq.'Y'.or.answer.eq.'y' )then
410   39   write(*,*) '  Enter the output filename : '
411        read(*,40) filout
412   40   format(a30)
413        open(3,file=filout,status='new',err=99)
414        goto 98
415c
416   99   write(*,*) '  File already exists, try again.'
417        go to 39
418c
419   98   continue
420        write(3,*) '  '
421        write(3,49) elips
422        finv=1.d0/f
423        b=a*(1.d0-f)
424        write(3,50) a,b,finv
425        write(*,*) '  Enter the First Station name : '
426        read(*,40) name1
427        write(*,*) '  Enter the Second Station name : '
428        read(*,40) name2
429c
430   41   format('  First  Station : ',a30)
431   42   format('  Second Station : ',a30)
432   84   format('  Error: First  Station ... Invalid Latitude  ')
433   85   format('  Error: First  Station ... Invalid Longitude ')
434   86   format('  Error: Second Station ... Invalid Latitude  ')
435   87   format('  Error: Second Station ... Invalid Longitude ')
436   88   format(1x,65(1h*))
437   91   format('         DD(0-89) MM(0-59) SS(0-59.999...)  ')
438   92   format('         DDD(0-359) MM(0-59) SS(0-59.999...)  ')
439c
440        write(3,*) '  '
441        write(3,41) name1
442        write(3,*) '  ---------------- '
443
444        if( irlat1.eq.0 )then
445          write(3,18) ld1, lm1, sl1, d_ns1
446        else
447          write(3,88)
448          write(3,84)
449          write(3,91)
450          write(3,88)
451        endif
452c
453        if( irlon1.eq.0 )then
454          write(3,19) lod1,lom1,slo1,d_ew1
455        else
456          write(3,88)
457          write(3,85)
458          write(3,92)
459          write(3,88)
460        endif
461c
462        write(3,*) '  '
463        write(3,42) name2
464        write(3,*) '  ---------------- '
465c
466        if( irlat2.eq.0 )then
467          write(3,18) ld2, lm2, sl2, d_ns2
468        else
469          write(3,88)
470          write(3,86)
471          write(3,91)
472          write(3,88)
473        endif
474c
475        if( irlon2.eq.0 )then
476          write(3,19) lod2,lom2,slo2,d_ew2
477        else
478          write(3,88)
479          write(3,87)
480          write(3,92)
481          write(3,88)
482        endif
483c
484        write(3,*) '  '
485        write(3,34) iaz1,maz1,saz1
486        write(3,35) iaz2,maz2,saz2
487        write(3,32) edist,dmt
488        write(3,*) '  '
489      endif
490c
491   80 write(*,*) '  '
492      write(*,*) '  '
493      write(*,*) '  '
494      write(*,*) '  1) Another inverse, different ellipsoid.'
495      write(*,*) '  2) Same ellipsoid, two new stations.'
496      write(*,*) '  3) Same ellipsoid, same First Station.'
497      write(*,*) '  4) Done.'
498      write(*,*) '  '
499      write(*,*) '  Enter choice : '
500      read(*,10) answer
501c
502      if(     answer.eq.'1' )then
503        goto 1
504      elseif( answer.eq.'2' )then
505        goto 15
506      elseif( answer.eq.'3' )then
507        goto 20
508      else
509        write(*,*) '  Thats all, folks!'
510      endif
511
512c     stop
513      end
514