cb::inverse
c
program inverse
c
c********1*********2*********3*********4*********5*********6*********7**
c
c name: inverse
c version: 200208.19
c author: stephen j. frakes
c purpose: to compute a geodetic inverse
c and then display output information
c
c input parameters:
c -----------------
c
c output parameters:
c ------------------
c
c local variables and constants:
c ------------------------------
c answer user prompt response
c b semiminor axis polar (in meters)
c baz azimuth back (in radians)
c buff input char buffer array
c dd,dm,ds temporary values for degrees, minutes, seconds
c dlon temporary value for difference in longitude (radians)
c dmt,d_mt char constants for meter units
c edist ellipsoid distance (in meters)
c elips ellipsoid choice
c esq eccentricity squared for reference ellipsoid
c faz azimuth forward (in radians)
c filout output file name
c finv reciprocal flattening
c hem hemisphere flag for lat & lon entry
c ierror error condition flag with d,m,s conversion
c lgh length of buff() array
c option user prompt response
c r1,r2 temporary variables
c ss temporary variable
c tol tolerance for conversion of seconds
c
c name1 name of station one
c ld1,lm1,sl1 latitude sta one - degrees,minutes,seconds
c ald1,alm1,sl1 latitude sta one - degrees,minutes,seconds
c lat1sn latitude sta one - sign (+/- 1)
c d_ns1 latitude sta one - char ('N','S')
c lod1,lom1,slo1 longitude sta one - degrees,minutes,seconds
c alod1,alom1,slo1 longitude sta one - degrees,minutes,seconds
c lon1sn longitude sta one - sign (+/- 1)
c d_ew1 longitude sta one - char ('E','W')
c iaz1,maz1,saz1 forward azimuth - degrees,minutes,seconds
c isign1 forward azimuth - flag (+/- 1)
c glat1,glon1 station one - (lat & lon in radians )
c p1,e1 standpoint one - (lat & lon in radians )
c
c name2 name of station two
c ld2,lm2,sl2 latitude sta two - degrees,minutes,seconds
c ald2,alm2,sl2 latitude sta two - degrees,minutes,seconds
c lat2sn latitude sta two - sign (+/- 1)
c d_ns2 latitude sta two - char ('N','S')
c lod2,lom2,slo2 longitude sta two - degrees,minutes,seconds
c alod2,alom2,slo2 longitude sta one - degrees,minutes,seconds
c lon2sn longitude sta two - sign (+/- 1)
c d_ew2 longitude sta two - char ('E','W')
c iaz2,maz2,saz2 back azimuth - degrees,minutes,seconds
c isign2 back azimuth - flag (+/- 1)
c glat2,glon2 station two - (lat & lon in radians )
c p2,e2 forepoint two - (lat & lon in radians )
c
c global variables and constants:
c -------------------------------
c a semimajor axis equatorial (in meters)
c f flattening
c pi constant 3.14159....
c rad constant 180.0/pi
c
c this module called by: n/a
c
c this module calls: elipss, getrad, inver1, todmsp
c gethem, trim, bufdms, gvalr8, gvali4, fixdms, gpnhri
c datan, write, read, dabs, open, stop
c
c include files used: n/a
c
c common blocks used: const, elipsoid
c
c references: see comments within subroutines
c
c comments:
c
c********1*********2*********3*********4*********5*********6*********7**
c::modification history
c::1990mm.dd, sjf, creation of program
c::199412.15, bmt, creation of program on viper
c::200203.08, crs, modified by c.schwarz to correct spelling of Clarke
c::200207.15, rws, modified i/o & standardized program documentation
c:: added subs trim, bufdms, gethem, gvali4, gvalr8
c::200207.23, rws, replaced sub inver1 with gpnarc, gpnloa, gpnhri
c::200208.15, rws, fixed an error in bufdms
c:: - renamed ellips to elipss "common error" with dirct2
c:: - added FAZ & BAZ to printed output
c::200208.19, rws, added more error flags for web interface code
c:: - added logical nowebb
c::200208.xx, sjf, program version number 2.0
c********1*********2*********3*********4*********5*********6*********7**
ce::inverse
c
implicit double precision (a-h, o-z)
implicit integer (i-n)
c
logical nowebb
c
character*1 answer,option,dmt,buff(50),hem
character*6 d_ns1, d_ew1, d_ns2, d_ew2, d_mt
character*30 filout,name1,name2,elips
c
integer*4 ierror
integer*4 lgh
c
common/const/pi,rad
common/elipsoid/a,f
c
c ms_unix 0 = web version
c 1 = ms_dos or unix version
c
parameter ( ms_unix = 0 )
c
pi = 4.d0*datan(1.d0)
rad = 180.d0/pi
dmt = 'm'
d_mt = 'Meters'
c
if( ms_unix.eq.1 )then
nowebb = .true.
else
nowebb = .false.
endif
c
1 do 2 i=1,25
write(*,*) ' '
2 continue
c
5 write(*,*) ' Program Inverse - Version 2.0 '
write(*,*) ' '
write(*,*) ' Ellipsoid options : '
write(*,*) ' '
write(*,*) ' 1) GRS80 / WGS84 (NAD83) '
write(*,*) ' 2) Clarke 1866 (NAD27) '
write(*,*) ' 3) Any other ellipsoid '
write(*,*) ' '
write(*,*) ' Enter choice : '
read(*,10) option
10 format(a1)
c
if(option.eq.'1') then
a=6378137.d0
f=1.d0/298.25722210088d0
elips='GRS80 / WGS84 (NAD83)'
elseif(option.eq.'2') then
a=6378206.4d0
f=1.d0/294.9786982138d0
elips='Clarke 1866 (NAD27)'
elseif(option.eq.'3') then
call elipss (elips)
else
write(*,*) ' Enter 1, 2, or 3 ! Try again --'
goto 5
endif
c
esq = f*(2.0d0-f)
c
write(*,*) ' '
write(*,*) ' '
write(*,*) ' '
write(*,*) ' '
c
15 write(*,*) ' Enter First Station '
write(*,16)
16 format(18x,'(Separate D,M,S by blanks or commas)')
write(*,*) 'hDD MM SS.sssss Latitude : (h default = N )'
11 format(50a1)
c
22 hem='N'
read(*,11) buff
call trim (buff,lgh,hem)
call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
c
if( ierror.eq.0 )then
irlat1 = 0
else
irlat1 = 1
write(*,*) ' Invalid Latitude ... Please re-enter it '
write(*,*) ' '
if( nowebb )then
goto 22
endif
endif
c
ald1 = dd
alm1 = dm
sl1 = ds
c
if( hem.eq.'N' ) then
lat1sn = +1
else
lat1sn = -1
endif
c
write(*,*) 'hDDD MM SS.sssss Longitude : (h default = W )'
c
23 hem='W'
read(*,11) buff
call trim (buff,lgh,hem)
call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
c
if( ierror.eq.0 )then
irlon1 = 0
else
irlon1 = 1
write(*,*) ' Invalid Longitude ... Please re-enter it '
write(*,*) ' '
if( nowebb )then
goto 23
endif
endif
c
alod1 = dd
alom1 = dm
slo1 = ds
c
if( hem.eq.'E' ) then
lon1sn = +1
else
lon1sn = -1
endif
c
call getrad(ald1, alm1, sl1, lat1sn, glat1)
call getrad(alod1,alom1,slo1,lon1sn, glon1)
c
write(*,*) ' '
write(*,*) ' '
c
20 write(*,*) ' Enter Second Station '
write(*,16)
write(*,*) 'hDD MM SS.sssss Latitude : (h default = N )'
c
24 hem='N'
read(*,11) buff
call trim (buff,lgh,hem)
call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
c
if( ierror.eq.0 )then
irlat2 = 0
else
irlat2 = 1
write(*,*) ' Invalid Latitude ... Please re-enter it '
write(*,*) ' '
if( nowebb )then
goto 24
endif
endif
c
ald2 = dd
alm2 = dm
sl2 = ds
c
if( hem.eq.'N' ) then
lat2sn = +1
else
lat2sn = -1
endif
c
write(*,*) 'hDDD MM SS.sssss Longitude : (h default = W )'
c
25 hem='W'
read(*,11) buff
call trim (buff,lgh,hem)
call bufdms (buff,lgh,hem,dd,dm,ds,ierror)
c
if( ierror.eq.0 )then
irlon2 = 0
else
irlon2 = 1
write(*,*) ' Invalid Longitude ... Please re-enter it '
write(*,*) ' '
if( nowebb )then
goto 25
endif
endif
c
alod2 = dd
alom2 = dm
slo2 = ds
c
if( hem.eq.'E' )then
lon2sn = +1
else
lon2sn = -1
endif
c
call getrad(ald2, alm2, sl2, lat2sn, glat2)
call getrad(alod2,alom2,slo2,lon2sn, glon2)
c
p1 = glat1
e1 = glon1
p2 = glat2
e2 = glon2
c
if( e1.lt.0.0d0 )then
e1 = e1+2.0d0*pi
endif
if( e2.lt.0.0d0 )then
e2 = e2+2.0d0*pi
endif
c
c compute the geodetic inverse
c
c ************************************************************
c * replaced subroutine inver1 with gpnhri
c *
c * call inver1 (glat1,glon1,glat2,glon2,faz,baz,edist)
c *
c ************************************************************
c
call gpnhri (a,f,esq,pi,p1,e1,p2,e2,faz,baz,edist)
c
c check for a non distance ... p1,e1 & p2,e2 equal zero ?
c
if( edist.lt.0.00005d0 )then
faz = 0.0d0
baz = 0.0d0
endif
c
c set the tolerance (in seconds) for the azimuth conversion
c
tol = 0.00005d0
c
call todmsp(faz,iaz1,maz1,saz1,isign1)
if(isign1.lt.0) then
iaz1=359-iaz1
maz1=59-maz1
saz1=60.d0-saz1
endif
call fixdms ( iaz1, maz1, saz1, tol )
c
call todmsp(baz,iaz2,maz2,saz2,isign2)
if(isign2.lt.0) then
iaz2=359-iaz2
maz2=59-maz2
saz2=60.d0-saz2
endif
call fixdms ( iaz2, maz2, saz2, tol )
c
call todmsp(glat1, ld1, lm1, sl1, lat1sn)
call todmsp(glon1, lod1, lom1, slo1, lon1sn)
call todmsp(glat2, ld2, lm2, sl2, lat2sn)
call todmsp(glon2, lod2, lom2, slo2, lon2sn)
c
call hem_ns ( lat1sn, d_ns1 )
call hem_ew ( lon1sn, d_ew1 )
call hem_ns ( lat2sn, d_ns2 )
call hem_ew ( lon2sn, d_ew2 )
c
write(*,*) ' '
write(*,*) ' '
write(*,*) ' '
write(*,*) ' '
write(*,*) ' '
write(*,49) elips
49 format(' Ellipsoid : ',a30)
finv=1.d0/f
b=a*(1.d0-f)
write(*,50) a,b,finv
50 format(' Equatorial axis, a = ',f15.4,/,
* ' Polar axis, b = ',f15.4,/,
* ' Inverse flattening, 1/f = ',f16.11)
c
18 format(' LAT = ',i3,1x,i2,1x,f8.5,1x,a6)
19 format(' LON = ',i3,1x,i2,1x,f8.5,1x,a6)
c
write(*,*) ' '
write(*,*) ' First Station : '
write(*,*) ' ---------------- '
write(*,18) ld1, lm1, sl1, d_ns1
write(*,19) lod1,lom1,slo1,d_ew1
c
write(*,*) ' '
write(*,*) ' Second Station : '
write(*,*) ' ---------------- '
write(*,18) ld2, lm2, sl2, d_ns2
write(*,19) lod2,lom2,slo2,d_ew2
c
32 format(' Ellipsoidal distance S = ',f14.4,1x,a1)
34 format(' Forward azimuth FAZ = ',i3,1x,i2,1x,f7.4,
1 ' From North')
35 format(' Back azimuth BAZ = ',i3,1x,i2,1x,f7.4,
1 ' From North')
c
write(*,*) ' '
write(*,34) iaz1,maz1,saz1
write(*,35) iaz2,maz2,saz2
write(*,32) edist,dmt
write(*,*) ' '
write(*,*) ' Do you want to save this output into a file (y/n)?'
read(*,10) answer
c
if( answer.eq.'Y'.or.answer.eq.'y' )then
39 write(*,*) ' Enter the output filename : '
read(*,40) filout
40 format(a30)
open(3,file=filout,status='new',err=99)
goto 98
c
99 write(*,*) ' File already exists, try again.'
go to 39
c
98 continue
write(3,*) ' '
write(3,49) elips
finv=1.d0/f
b=a*(1.d0-f)
write(3,50) a,b,finv
write(*,*) ' Enter the First Station name : '
read(*,40) name1
write(*,*) ' Enter the Second Station name : '
read(*,40) name2
c
41 format(' First Station : ',a30)
42 format(' Second Station : ',a30)
84 format(' Error: First Station ... Invalid Latitude ')
85 format(' Error: First Station ... Invalid Longitude ')
86 format(' Error: Second Station ... Invalid Latitude ')
87 format(' Error: Second Station ... Invalid Longitude ')
88 format(1x,65(1h*))
91 format(' DD(0-89) MM(0-59) SS(0-59.999...) ')
92 format(' DDD(0-359) MM(0-59) SS(0-59.999...) ')
c
write(3,*) ' '
write(3,41) name1
write(3,*) ' ---------------- '
if( irlat1.eq.0 )then
write(3,18) ld1, lm1, sl1, d_ns1
else
write(3,88)
write(3,84)
write(3,91)
write(3,88)
endif
c
if( irlon1.eq.0 )then
write(3,19) lod1,lom1,slo1,d_ew1
else
write(3,88)
write(3,85)
write(3,92)
write(3,88)
endif
c
write(3,*) ' '
write(3,42) name2
write(3,*) ' ---------------- '
c
if( irlat2.eq.0 )then
write(3,18) ld2, lm2, sl2, d_ns2
else
write(3,88)
write(3,86)
write(3,91)
write(3,88)
endif
c
if( irlon2.eq.0 )then
write(3,19) lod2,lom2,slo2,d_ew2
else
write(3,88)
write(3,87)
write(3,92)
write(3,88)
endif
c
write(3,*) ' '
write(3,34) iaz1,maz1,saz1
write(3,35) iaz2,maz2,saz2
write(3,32) edist,dmt
write(3,*) ' '
endif
c
80 write(*,*) ' '
write(*,*) ' '
write(*,*) ' '
write(*,*) ' 1) Another inverse, different ellipsoid.'
write(*,*) ' 2) Same ellipsoid, two new stations.'
write(*,*) ' 3) Same ellipsoid, same First Station.'
write(*,*) ' 4) Done.'
write(*,*) ' '
write(*,*) ' Enter choice : '
read(*,10) answer
c
if( answer.eq.'1' )then
goto 1
elseif( answer.eq.'2' )then
goto 15
elseif( answer.eq.'3' )then
goto 20
else
write(*,*) ' Thats all, folks!'
endif
c stop
end
subroutine bufdms (buff,lgh,hem,dd,dm,ds,ierror)
implicit double precision (a-h, o-z)
implicit integer (i-n)
c
logical done,flag
c
character*1 buff(*),abuf(21)
character*1 ch
character*1 hem
integer*4 ll,lgh
integer*4 i4,id,im,is,icond,ierror
real*8 x(5)
c
c set the "error flag"
c
ierror = 0
icond = 0
c
c set defaults for dd,dm,ds
c
dd = 0.0d0
dm = 0.0d0
ds = 0.0d0
c
c set default limits for "hem" flag
c
if( hem.eq.'N' .or. hem.eq.'S' )then
ddmax = 90.0d0
elseif( hem.eq.'E' .or. hem.eq.'W' )then
ddmax = 360.0d0
elseif( hem.eq.'A' )then
ddmax = 360.0d0
elseif( hem.eq.'Z' )then
ddmax = 180.0d0
elseif( hem.eq.'*' )then
ddmax = 0.0d0
ierror = 1
else
ddmax = 360.0d0
endif
c
do 1 i=1,5
x(i) = 0.0d0
1 continue
c
icolon = 0
ipoint = 0
icount = 0
flag = .true.
jlgh = lgh
c
do 2 i=1,jlgh
if( buff(i).eq.':' )then
icolon = icolon+1
endif
if( buff(i).eq.'.' )then
ipoint = ipoint+1
flag = .false.
endif
if( flag )then
icount = icount+1
endif
2 continue
c
if( ipoint.eq.1 .and. icolon.eq.0 )then
c
c load temp buffer
c
do 3 i=1,jlgh
abuf(i) = buff(i)
3 continue
abuf(jlgh+1) = '$'
ll = jlgh
c
call gvalr8 (abuf,ll,r8,icond)
c
if( icount.ge.5 )then
c
c value is a packed decimal of ==> DDMMSS.sssss
c
ss = r8/10000.0d0
id = idint( ss )
c
r8 = r8-10000.0d0*dble(float(id))
ss = r8/100.0d0
im = idint( ss )
c
r8 = r8-100.0d0*dble(float(im))
else
c
c value is a decimal of ==> .xx X.xxx X.
c
id = idint( r8 )
r8 = (r8-id)*60.0d0
im = idint( r8 )
r8 = (r8-im)*60.0d0
endif
c
c account for rounding error
c
is = idnint( r8*1.0d5 )
if( is.ge.6000000 )then
r8 = 0.0d0
im = im+1
endif
c
if( im.ge.60 )then
im = 0
id = id+1
endif
c
dd = dble( float( id ) )
dm = dble( float( im ) )
ds = r8
else
c
c buff() value is a d,m,s of ==> NN:NN:XX.xxx
c
k = 0
next = 1
done = .false.
ie = jlgh
c
do 100 j=1,5
ib = next
do 90 i=ib,ie
ch = buff(i)
last = i
if( i.eq.jlgh .or. ch.eq.':' )then
if( i.eq.jlgh )then
done = .true.
endif
if( ch.eq.':' )then
last = i-1
endif
goto 91
endif
90 continue
goto 98
c
91 ipoint = 0
ik = 0
do 92 i=next,last
ik = ik+1
ch = buff(i)
if( ch.eq.'.' )then
ipoint = ipoint+1
endif
abuf(ik) = buff(i)
92 continue
abuf(ik+1) = '$'
c
ll = ik
if( ipoint.eq.0 )then
call gvali4 (abuf,ll,i4,icond)
r8 = dble(float( i4 ))
else
call gvalr8 (abuf,ll,r8,icond)
endif
c
k = k+1
x(k) = r8
c
98 if( done )then
goto 101
endif
c
next = last
99 next = next+1
if( buff(next).eq.':' )then
goto 99
endif
100 continue
c
c load dd,dm,ds
c
101 if( k.ge.1 )then
dd = x(1)
endif
c
if( k.ge.2 )then
dm = x(2)
endif
c
if( k.ge.3 )then
ds = x(3)
endif
endif
c
if( dd.gt.ddmax .or.
1 dm.ge.60.0d0 .or.
1 ds.ge.60.0d0 )then
ierror = 1
dd = 0.0d0
dm = 0.0d0
ds = 0.0d0
endif
c
if( icond.ne.0 )then
ierror = 1
endif
c
return
end
subroutine elipss (elips)
implicit double precision(a-h,o-z)
character*1 answer
character*30 elips
common/elipsoid/a,f
write(*,*) ' Other Ellipsoids.'
write(*,*) ' -----------------'
write(*,*) ' '
write(*,*) ' A) Airy 1858'
write(*,*) ' B) Airy Modified'
write(*,*) ' C) Australian National'
write(*,*) ' D) Bessel 1841'
write(*,*) ' E) Clarke 1880'
write(*,*) ' F) Everest 1830'
write(*,*) ' G) Everest Modified'
write(*,*) ' H) Fisher 1960'
write(*,*) ' I) Fisher 1968'
write(*,*) ' J) Hough 1956'
write(*,*) ' K) International (Hayford)'
write(*,*) ' L) Krassovsky 1938'
write(*,*) ' M) NWL-9D (WGS 66)'
write(*,*) ' N) South American 1969'
write(*,*) ' O) Soviet Geod. System 1985'
write(*,*) ' P) WGS 72'
write(*,*) ' Q-Z) User defined.'
write(*,*) ' '
write(*,*) ' Enter choice : '
read(*,10) answer
10 format(a1)
c
if(answer.eq.'A'.or.answer.eq.'a') then
a=6377563.396d0
f=1.d0/299.3249646d0
elips='Airy 1858'
elseif(answer.eq.'B'.or.answer.eq.'b') then
a=6377340.189d0
f=1.d0/299.3249646d0
elips='Airy Modified'
elseif(answer.eq.'C'.or.answer.eq.'c') then
a=6378160.d0
f=1.d0/298.25d0
elips='Australian National'
elseif(answer.eq.'D'.or.answer.eq.'d') then
a=6377397.155d0
f=1.d0/299.1528128d0
elips='Bessel 1841'
elseif(answer.eq.'E'.or.answer.eq.'e') then
a=6378249.145d0
f=1.d0/293.465d0
elips='Clarke 1880'
elseif(answer.eq.'F'.or.answer.eq.'f') then
a=6377276.345d0
f=1.d0/300.8017d0
elips='Everest 1830'
elseif(answer.eq.'G'.or.answer.eq.'g') then
a=6377304.063d0
f=1.d0/300.8017d0
elips='Everest Modified'
elseif(answer.eq.'H'.or.answer.eq.'h') then
a=6378166.d0
f=1.d0/298.3d0
elips='Fisher 1960'
elseif(answer.eq.'I'.or.answer.eq.'i') then
a=6378150.d0
f=1.d0/298.3d0
elips='Fisher 1968'
elseif(answer.eq.'J'.or.answer.eq.'j') then
a=6378270.d0
f=1.d0/297.d0
elips='Hough 1956'
elseif(answer.eq.'K'.or.answer.eq.'k') then
a=6378388.d0
f=1.d0/297.d0
elips='International (Hayford)'
elseif(answer.eq.'L'.or.answer.eq.'l') then
a=6378245.d0
f=1.d0/298.3d0
elips='Krassovsky 1938'
elseif(answer.eq.'M'.or.answer.eq.'m') then
a=6378145.d0
f=1.d0/298.25d0
elips='NWL-9D (WGS 66)'
elseif(answer.eq.'N'.or.answer.eq.'n') then
a=6378160.d0
f=1.d0/298.25d0
elips='South American 1969'
elseif(answer.eq.'O'.or.answer.eq.'o') then
a=6378136.d0
f=1.d0/298.257d0
elips='Soviet Geod. System 1985'
elseif(answer.eq.'P'.or.answer.eq.'p') then
a=6378135.d0
f=1.d0/298.26d0
elips='WGS 72'
else
elips = 'User defined.'
c
write(*,*) ' Enter Equatorial axis, a : '
read(*,*) a
a = dabs(a)
c
write(*,*) ' Enter either Polar axis, b or '
write(*,*) ' Reciprocal flattening, 1/f : '
read(*,*) ss
ss = dabs(ss)
c
f = 0.0d0
if( 200.0d0.le.ss .and. ss.le.310.0d0 )then
f = 1.d0/ss
elseif( 6000000.0d0.lt.ss .and. ss.lt.a )then
f = (a-ss)/a
else
elips = 'Error: default GRS80 used.'
a = 6378137.0d0
f = 1.0d0/298.25722210088d0
endif
endif
c
return
end
subroutine fixdms (ideg, min, sec, tol )
c
implicit double precision (a-h, o-z)
implicit integer (i-n)
c
c test for seconds near 60.0-tol
c
if( sec.ge.( 60.0d0-tol ) )then
sec = 0.0d0
min = min+1
endif
c
c test for minutes near 60
c
if( min.ge.60 )then
min = 0
ideg = ideg+1
endif
c
c test for degrees near 360
c
if( ideg.ge.360 )then
ideg = 0
endif
c
return
end
subroutine hem_ns ( lat_sn, hem )
implicit integer (i-n)
character*6 hem
c
if( lat_sn.eq.1 ) then
hem = 'North '
else
hem = 'South '
endif
c
return
end
subroutine hem_ew ( lon_sn, hem )
implicit integer (i-n)
character*6 hem
c
if( lon_sn.eq.1 ) then
hem = 'East '
else
hem = 'West '
endif
c
return
end
subroutine getrad(d,m,sec,isign,val)
*** comvert deg, min, sec to radians
implicit double precision(a-h,j-z)
common/const/pi,rad
val=(d+m/60.d0+sec/3600.d0)/rad
val=dble(isign)*val
return
end
CB::GPNARC
C
SUBROUTINE GPNARC (AMAX,FLAT,ESQ,PI,P1,P2,ARC)
C
C********1*********2*********3*********4*********5*********6*********7*
C
C NAME: GPNARC
C VERSION: 200005.26
C WRITTEN BY: ROBERT (Sid) SAFFORD
C PURPOSE: SUBROUTINE TO COMPUTE THE LENGTH OF A MERIDIONAL ARC
C BETWEEN TWO LATITUDES
C
C INPUT PARAMETERS:
C -----------------
C AMAX SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID
C FLAT FLATTENING (0.0033528 ... )
C ESQ ECCENTRICITY SQUARED FOR REFERENCE ELLIPSOID
C PI 3.14159...
C P1 LAT STATION 1
C P2 LAT STATION 2
C
C OUTPUT PARAMETERS:
C ------------------
C ARC GEODETIC DISTANCE
C
C LOCAL VARIABLES AND CONSTANTS:
C ------------------------------
C GLOBAL VARIABLES AND CONSTANTS:
C -------------------------------
C
C MODULE CALLED BY: GENERAL
C
C THIS MODULE CALLS:
C LLIBFORE/ OPEN, CLOSE, READ, WRITE, INQUIRE
C DABS, DBLE, FLOAT, IABS, CHAR, ICHAR
C
C INCLUDE FILES USED:
C COMMON BLOCKS USED:
C
C REFERENCES: Microsoft FORTRAN 4.10 Optimizing Compiler, 1988
C MS-DOS Operating System
C COMMENTS:
C********1*********2*********3*********4*********5*********6*********7*
C::MODIFICATION HISTORY
C::197507.05, RWS, VER 00 TENCOL RELEASED FOR FIELD USE
C::198311.20, RWS, VER 01 MTEN RELEASED TO FIELD
C::198411.26, RWS, VER 07 MTEN2 RELEASED TO FIELD
C::1985xx.xx, RWS, CODE CREATED
C::198506.10, RWS, WRK ENHANCEMENTS RELEASED TO FIELD
C::198509.01, RWS, VER 11 MTEN3 RELEASED TO FIELD
C::198512.18, RWS, CODE MODIFIED FOR MTEN3
C::198708.10, RWS, CODE MODIFIED TO USE NEW MTEN4 GPN RECORD FORMAT
C::199112.31, RWS, VER 20 MTEN4 RELEASED TO FIELD
C::200001.13, RWS, VER 21 MTEN4 RELEASED TO FIELD
C::200005.26, RWS, CODE RESTRUCTURED & DOCUMENTATION ADDED
C::200012.31, RWS, VER 23 MTEN5 RELEASED
C********1*********2*********3*********4*********5*********6*********7*
CE::GPNARC
C ---------------------------
C M T E N (VERSION 3)
C M T E N (VERSION 5.23)
C ---------------------------
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LOGICAL FLAG
C
DATA TT/5.0D-15/
C
C CHECK FOR A 90 DEGREE LOOKUP
C
FLAG = .FALSE.
C
S1 = DABS(P1)
S2 = DABS(P2)
C
IF( (PI/2.0D0-TT).LT.S2 .AND. S2.LT.(PI/2.0D0+TT) )THEN
FLAG = .TRUE.
ENDIF
C
IF( S1.GT.TT )THEN
FLAG = .FALSE.
ENDIF
C
DA = (P2-P1)
S1 = 0.0D0
S2 = 0.0D0
C
C COMPUTE THE LENGTH OF A MERIDIONAL ARC BETWEEN TWO LATITUDES
C
E2 = ESQ
E4 = E2*E2
E6 = E4*E2
E8 = E6*E2
EX = E8*E2
C
T1 = E2*(003.0D0/4.0D0)
T2 = E4*(015.0D0/64.0D0)
T3 = E6*(035.0D0/512.0D0)
T4 = E8*(315.0D0/16384.0D0)
T5 = EX*(693.0D0/131072.0D0)
C
A = 1.0D0+T1+3.0D0*T2+10.0D0*T3+35.0D0*T4+126.0D0*T5
C
IF( FLAG )THEN
GOTO 1
ENDIF
C
B = T1+4.0D0*T2+15.0D0*T3+56.0D0*T4+210.0D0*T5
C = T2+06.0D0*T3+28.0D0*T4+120.0D0*T5
D = T3+08.0D0*T4+045.0D0*T5
E = T4+010.0D0*T5
F = T5
C
DB = DSIN(P2*2.0D0)-DSIN(P1*2.0D0)
DC = DSIN(P2*4.0D0)-DSIN(P1*4.0D0)
DD = DSIN(P2*6.0D0)-DSIN(P1*6.0D0)
DE = DSIN(P2*8.0D0)-DSIN(P1*8.0D0)
DF = DSIN(P2*10.0D0)-DSIN(P1*10.0D0)
C
C COMPUTE THE S2 PART OF THE SERIES EXPANSION
C
S2 = -DB*B/2.0D0+DC*C/4.0D0-DD*D/6.0D0+DE*E/8.0D0-DF*F/10.0D0
C
C COMPUTE THE S1 PART OF THE SERIES EXPANSION
C
1 S1 = DA*A
C
C COMPUTE THE ARC LENGTH
C
ARC = AMAX*(1.0D0-ESQ)*(S1+S2)
C
RETURN
END
cb::gpnhri
c
subroutine gpnhri (a,f,esq,pi,p1,e1,p2,e2,az1,az2,s)
c
c********1*********2*********3*********4*********5*********6*********7*
c
c name: gpnhri
c version: 200208.09
c written by: robert (sid) safford
c purpose: subroutine to compute helmert rainsford inverse problem
c
c solution of the geodetic inverse problem after t. vincenty
c modified rainsford's method with helmert's elliptical terms
c effective in any azimuth and at any distance short of antipocal
c from/to stations must not be the geographic pole.
c parameter a is the semi-major axis of the reference ellipsoid
c finv=1/f is the inverse flattening of the reference ellipsoid
c latitudes and longitudes in radians positive north and west
c forward and back azimuths returned in radians clockwise from south
c geodesic distance s returned in units of semi-major axis a
c programmed for ibm 360-195 09/23/75
c
c note - note - note -
c 1. do not use for meridional arcs and be careful on the equator.
c 2. azimuths are from north(+) clockwise and
c 3. longitudes are positive east(+)
c
c input parameters:
c -----------------
c a semi-major axis of reference ellipsoid meters
c f flattening (0.0033528...)
c esq eccentricity squared
c pi 3.14159...
c p1 lat station 1 radians
c e1 lon station 1 radians
c p2 lat station 2 radians
c e2 lon station 2 radians
c
c output parameters:
c ------------------
c az1 azi at sta 1 -> sta 2 radians
c az2 azi at sta 2 -> sta 1 radians
c s geodetic dist between sta(s) 1 & 2 meters
c
c local variables and constants:
c ------------------------------
c aa constant from subroutine gpnloa
c alimit equatorial arc distance along the equator (radians)
c arc meridional arc distance latitude p1 to p2 (in meters)
c az1 azimuth forward (in radians)
c az2 azimuth back (in radians)
c bb constant from subroutine gpnloa
c dlon temporary value for difference in longitude (radians)
c equ equatorial distance (in meters)
c r1,r2 temporary variables
c s ellipsoid distance (in meters)
c sms equatorial - geodesic distance (S - s) "Sms"
c ss temporary variable
c tol0 tolerance for checking computation value
c tol1 tolerance for checking a real zero value
c tol2 tolerance for close to zero value
c twopi two times constant pi
c
c global variables and constants:
c -------------------------------
c
c module called by: general
c
c this module calls: gpnarc, gpnloa
c llibfore/ dsin, dcos, dsqrt, dabs, datan2, write
c
c include files used:
c common blocks used:
c
c references: microsoft fortran 4.10 optimizing compiler, 1988
c ms-dos operating system
c comments:
c********1*********2*********3*********4*********5*********6*********7*
c::modification history
c::197507.05, rws, ver 00 tencol released for field use
c::198311.20, rws, ver 01 mten released to field
c::198411.26, rws, ver 07 mten2 released to field
c::198506.10, rws, wrk enhancements released to field
c::198507.22, rws, code modified for mten3
c::198509.01, rws, ver 11 mten3 released to field
c::198708.10, rws, code modified to use new mten4 gpn record format
c::199112.31, rws, ver 20 mten4 released to field
c::200001.13, rws, ver 21 mten4 released to field
c::200005.26, rws, code restructured & documentation added
c::200012.31, rws, ver 23 mten5 released
c::200104.09, rws, code added to calblin program
c::200208.09, rws, code added subroutines gpnarc & gpnloa
c********1*********2*********3*********4*********5*********6*********7*
ce::gpnhri
c -------------------------------
c m t e n (version 3)
c (version 4.22)
c (version 5.23)
c -------------------------------
c
implicit real*8 (a-h,o-z)
c
data tol0 /5.0d-15/
data tol1 /5.0d-14/
data tol2 /7.0d-03/
c
twopi = 2.0d0*pi
c
c test the longitude difference with tol1
c tol1 is approximately 0.000000001 arc seconds
c
ss = e2-e1
if( dabs(ss).lt.tol1 )then
e2 = e2+tol1
write(*,*) ' longitudal difference is near zero '
c
r2 = p2
r1 = p1
call gpnarc ( a, f, esq, pi, r1, r2, arc )
s = dabs( arc )
c
if( p2.gt.p1 )then
az1 = 0.0d0
az2 = pi
else
az1 = pi
az2 = 0.0d0
endif
return
endif
c
c test for longitude over 180 degrees
c
dlon = e2-e1
c
if( dlon.ge.0.0d0 )then
if( pi.le.dlon .and. dlon.lt.twopi )then
dlon = dlon-twopi
endif
else
ss = dabs(dlon)
if( pi.le.ss .and. ss.lt.twopi )then
dlon = dlon+twopi
endif
endif
c
ss = dabs( dlon )
if( ss.gt.pi )then
c:: write(*,*) ' '
c:: write(*,*) ' Longitude difference over 180 degrees '
c:: write(*,*) ' Turn it around '
ss = twopi-ss
endif
c
c compute the limit in longitude (alimit), it is equal
c to twice the distance from the equator to the pole,
c as measured along the equator (east/ewst)
c
alimit = pi*(1.0d0-f)
c
c test for anti-nodal difference
c
if( ss.ge.alimit )then
r1 = dabs(p1)
r2 = dabs(p2)
c
c latitudes r1 & r2 are not near the equator
c
if( r1.gt.tol2 .and. r2.gt.tol2 )then
goto 60
endif
c
c longitude difference is greater than lift-off point
c now check to see if "both" r1 & r2 are on equator
c
if( r1.lt.tol1 .and. r2.gt.tol2 )then
goto 60
endif
if( r2.lt.tol1 .and. r1.gt.tol2 )then
goto 60
endif
c
c check for either r1 or r2 just off the equator but < tol2
c
if( r1.gt.tol1. or. r2.gt.tol1 )then
az1 = 0.0d0
az2 = 0.0d0
s = 0.0d0
return
endif
c
c compute the azimuth to anti-nodal point
c
c:: write(*,*) ' '
c:: write(*,*) ' Longitude difference beyond lift-off point '
c:: write(*,*) ' '
c
call gpnloa (a,f,esq,pi,dlon,az1,az2,aa,bb,sms)
c
c compute the equatorial distance & geodetic
c
equ = a*dabs(dlon)
s = equ-sms
return
endif
c
60 continue
c
f0 = (1.0d0-f)
b = a*f0
epsq = esq/(1.0d0-esq)
f2 = f*f
f3 = f*f2
f4 = f*f3
c
c the longitude difference
c
dlon = e2-e1
ab = dlon
kount = 0
c
c the reduced latitudes
c
u1 = f0*dsin(p1)/dcos(p1)
u2 = f0*dsin(p2)/dcos(p2)
c
u1 = datan(u1)
u2 = datan(u2)
c
su1 = dsin(u1)
cu1 = dcos(u1)
c
su2 = dsin(u2)
cu2 = dcos(u2)
c
c counter for the iteration operation
c
1 kount = kount+1
c
clon = dcos(ab)
slon = dsin(ab)
c
csig = su1*su2+cu1*cu2*clon
ssig = dsqrt((slon*cu2)**2+(su2*cu1-su1*cu2*clon)**2)
c
sig = datan2(ssig,csig)
sinalf=cu1*cu2*slon/ssig
c
w = (1.0d0-sinalf*sinalf)
t4 = w*w
t6 = w*t4
c
c the coefficients of type a
c
ao = f-f2*(1.0d0+f+f2)*w/4.0d0+3.0d0*f3*(1.0d0+
1 9.0d0*f/4.0d0)*t4/16.0d0-25.0d0*f4*t6/128.0d0
a2 = f2*(1.0d0+f+f2)*w/4.0d0-f3*(1.0d0+9.0d0*f/4.0d0)*t4/4.0d0+
1 75.0d0*f4*t6/256.0d0
a4 = f3*(1.0d0+9.0d0*f/4.0d0)*t4/32.0d0-15.0d0*f4*t6/256.0d0
a6 = 5.0d0*f4*t6/768.0d0
c
c the multiple angle functions
c
qo = 0.0d0
if( w.gt.tol0 )then
qo = -2.0d0*su1*su2/w
endif
c
q2 = csig+qo
q4 = 2.0d0*q2*q2-1.0d0
q6 = q2*(4.0d0*q2*q2-3.0d0)
r2 = 2.0d0*ssig*csig
r3 = ssig*(3.0d0-4.0d0*ssig*ssig)
c
c the longitude difference
c
s = sinalf*(ao*sig+a2*ssig*q2+a4*r2*q4+a6*r3*q6)
xz = dlon+s
c
xy = dabs(xz-ab)
ab = dlon+s
c
if( xy.lt.0.5d-13 )then
goto 4
endif
c
if( kount.le.7 )then
goto 1
endif
c
c the coefficients of type b
c
4 z = epsq*w
c
bo = 1.0d0+z*(1.0d0/4.0d0+z*(-3.0d0/64.0d0+z*(5.0d0/256.0d0-
1 z*175.0d0/16384.0d0)))
b2 = z*(-1.0d0/4.0d0+z*(1.0d0/16.0d0+z*(-15.0d0/512.0d0+
1 z*35.0d0/2048.0d0)))
b4 = z*z*(-1.0d0/128.0d0+z*(3.0d0/512.0d0-z*35.0d0/8192.0d0))
b6 = z*z*z*(-1.0d0/1536.0d0+z*5.0d0/6144.0d0)
c
c the distance in meters
c
s = b*(bo*sig+b2*ssig*q2+b4*r2*q4+b6*r3*q6)
c
c first compute the az1 & az2 for along the equator
c
if( dlon.gt.pi )then
dlon = (dlon-2.0d0*pi)
endif
c
if( dabs(dlon).gt.pi )then
dlon = (dlon+2.0d0*pi)
endif
c
az1 = pi/2.0d0
if( dlon.lt.0.0d0 )then
az1 = 3.0d0*az1
endif
c
az2 = az1+pi
if( az2.gt.2.0d0*pi )then
az2 = az2-2.0d0*pi
endif
c
c now compute the az1 & az2 for latitudes not on the equator
c
if( .not.(dabs(su1).lt.tol0 .and. dabs(su2).lt.tol0) )then
tana1 = slon*cu2/(su2*cu1-clon*su1*cu2)
tana2 = slon*cu1/(su1*cu2-clon*su2*cu1)
sina1 = sinalf/cu1
sina2 = -sinalf/cu2
c
c azimuths from north,longitudes positive east
c
az1 = datan2(sina1,sina1/tana1)
az2 = pi-datan2(sina2,sina2/tana2)
endif
c
if( az1.lt.0.0d0 )then
az1 = az1+2.0d0*pi
endif
c
if( az2.lt.0.0d0 )then
az2 = az2+2.0d0*pi
endif
c
return
end
CB::GPNLOA
C
SUBROUTINE GPNLOA (AMAX,FLAT,ESQ,PI,DL,AZ1,AZ2,AO,BO,SMS)
C
C********1*********2*********3*********4*********5*********6*********7*
C
C NAME: GPNLOA
C VERSION: 200005.26
C WRITTEN BY: ROBERT (Sid) SAFFORD
C PURPOSE: SUBROUTINE TO COMPUTE THE LIFF-OFF-AZIMUTH CONSTANTS
C
C INPUT PARAMETERS:
C -----------------
C AMAX SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID
C FLAT FLATTENING (0.0033528 ... )
C ESQ ECCENTRICITY SQUARED FOR REFERENCE ELLIPSOID
C PI 3.14159...
C DL LON DIFFERENCE
C AZ1 AZI AT STA 1 -> STA 2
C
C OUTPUT PARAMETERS:
C ------------------
C AZ2 AZ2 AT STA 2 -> STA 1
C AO CONST
C BO CONST
C SMS DISTANCE ... EQUATORIAL - GEODESIC (S - s) "SMS"
C
C LOCAL VARIABLES AND CONSTANTS:
C ------------------------------
C GLOBAL VARIABLES AND CONSTANTS:
C -------------------------------
C
C MODULE CALLED BY: GENERAL
C
C THIS MODULE CALLS:
C LLIBFORE/ DSIN, DCOS, DABS, DASIN
C
C INCLUDE FILES USED:
C COMMON BLOCKS USED:
C
C REFERENCES: Microsoft FORTRAN 4.10 Optimizing Compiler, 1988
C MS-DOS Operating System
C COMMENTS:
C********1*********2*********3*********4*********5*********6*********7*
C::MODIFICATION HISTORY
C::1985xx.xx, RWS, CODE CREATED
C::198506.10, RWS, WRK ENHANCEMENTS RELEASED TO FIELD
C::198509.01, RWS, VER 11 MTEN3 RELEASED TO FIELD
C::198512.18, RWS, CODE MODIFIED FOR MTEN3
C::198708.10, RWS, CODE MODIFIED TO USE NEW MTEN4 GPN RECORD FORMAT
C::199112.31, RWS, VER 20 MTEN4 RELEASED TO FIELD
C::200001.13, RWS, VER 21 MTEN4 RELEASED TO FIELD
C::200005.26, RWS, CODE RESTRUCTURED & DOCUMENTATION ADDED
C::200012.31, RWS, VER 23 MTEN5 RELEASED
C********1*********2*********3*********4*********5*********6*********7*
CE::GPNLOA
C ---------------------------
C M T E N (VERSION 3)
C (VERSION 4.22)
C (VERSION 5.23)
C ---------------------------
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DATA TT/5.0D-13/
C
DLON = DABS(DL)
CONS = (PI-DLON)/(PI*FLAT)
F = FLAT
C
C COMPUTE AN APPROXIMATE AZ
C
AZ = DASIN(CONS)
C
T1 = 1.0D0
T2 = (-1.0D0/4.0D0)*F*(1.0D0+F+F*F)
T4 = 3.0D0/16.0D0*F*F*(1.0D0+(9.0D0/4.0D0)*F)
T6 = (-25.0D0/128.0D0)*F*F*F
C
ITER = 0
1 ITER = ITER+1
S = DCOS(AZ)
C2 = S*S
C
C COMPUTE NEW AO
C
AO = T1 + T2*C2 + T4*C2*C2 + T6*C2*C2*C2
CS = CONS/AO
S = DASIN(CS)
IF( DABS(S-AZ).LT.TT )THEN
GOTO 2
ENDIF
C
AZ = S
IF( ITER.LE.6 )THEN
GOTO 1
ENDIF
C
2 AZ1 = S
IF( DL.LT.0.0D0 )THEN
AZ1 = 2.0D0*PI-AZ1
ENDIF
C
AZ2 = 2.0D0*PI-AZ1
C
C EQUATORIAL - GEODESIC (S - s) "SMS"
C
ESQP = ESQ/(1.0D0-ESQ)
S = DCOS(AZ1)
C
U2 = ESQP*S*S
U4 = U2*U2
U6 = U4*U2
U8 = U6*U2
C
T1 = 1.0D0
T2 = (1.0D0/4.0D0)*U2
T4 = (-3.0D0/64.0D0)*U4
T6 = (5.0D0/256.0D0)*U6
T8 = (-175.0D0/16384.0D0)*U8
C
BO = T1 + T2 + T4 + T6 + T8
S = DSIN(AZ1)
SMS = AMAX*PI*(1.0D0 - FLAT*DABS(S)*AO - BO*(1.0D0-FLAT))
C
RETURN
END
subroutine gvali4 (buff,ll,vali4,icond)
implicit integer (i-n)
c
logical plus,sign,done,error
character*1 buff(*)
character*1 ch
c
c integer*2 i
c integer*2 l1
c
integer*4 ich,icond
integer*4 ll
integer*4 vali4
c
l1 = ll
vali4 = 0
icond = 0
plus = .true.
sign = .false.
done = .false.
error = .false.
c
i = 0
10 i = i+1
if( i.gt.l1 .or. done )then
go to 1000
else
ch = buff(i)
ich = ichar( buff(i) )
endif
c
if( ch.eq.'+' )then
c
c enter on plus sign
c
if( sign )then
goto 150
else
sign = .true.
goto 10
endif
elseif( ch.eq.'-' )then
c
c enter on minus sign
c
if( sign )then
goto 150
else
sign = .true.
plus = .false.
goto 10
endif
elseif( ch.ge.'0' .and. ch.le.'9' )then
goto 100
elseif( ch.eq.' ' )then
c
c enter on space -- ignore leading spaces
c
if( .not.sign )then
goto 10
else
buff(i) = '0'
ich = 48
goto 100
endif
elseif( ch.eq.':' )then
c
c enter on colon -- ignore
c
if( .not.sign )then
goto 10
else
goto 1000
endif
elseif( ch.eq.'$' )then
c
c enter on dollar "$"
c
done = .true.
goto 10
else
c
c something wrong
c
goto 150
endif
c
c enter on numeric
c
100 vali4 = 10*vali4+(ich-48)
sign = .true.
goto 10
c
c treat illegal character
c
150 buff(i) = '0'
vali4 = 0
icond = 1
c
1000 if( .not.plus )then
vali4 = -vali4
endif
c
return
end
subroutine gvalr8 (buff,ll,valr8,icond)
implicit integer (i-n)
c
logical plus,sign,dpoint,done
c
character*1 buff(*)
character*1 ch
c
c integer*2 i, ip
c integer*2 l1
c integer*2 nn, num, n48
c
integer*4 ich,icond
integer*4 ll
c
real*8 ten
real*8 valr8
real*8 zero
c
data zero,ten/0.0d0,10.0d0/
c
n48 = 48
l1 = ll
icond = 0
valr8 = zero
plus = .true.
sign = .false.
dpoint = .false.
done = .false.
c
c start loop thru buffer
c
i = 0
10 i = i+1
if( i.gt.l1 .or. done )then
go to 1000
else
ch = buff(i)
nn = ichar( ch )
ich = nn
endif
c
if( ch.eq.'+' )then
c
c enter on plus sign
c
if( sign )then
goto 150
else
sign = .true.
goto 10
endif
elseif( ch.eq.'-' )then
c
c enter on minus sign
c
if( sign )then
goto 150
else
sign = .true.
plus = .false.
goto 10
endif
elseif( ch.eq.'.' )then
c
c enter on decimal point
c
ip = 0
sign = .true.
dpoint = .true.
goto 10
elseif( ch.ge.'0' .and. ch.le.'9' )then
goto 100
elseif( ch.eq.' ' )then
c
c enter on space
c
if( .not.sign )then
goto 10
else
buff(i) = '0'
ich = 48
goto 100
endif
elseif( ch.eq.':' .or. ch.eq.'$' )then
c
c enter on colon or "$" sign
c
done = .true.
goto 10
else
c
c something wrong
c
goto 150
endif
c
c enter on numeric
c
100 sign = .true.
if( dpoint )then
ip = ip+1
endif
c
num = ich
valr8 = ten*valr8+dble(float( num-n48 ))
goto 10
c
c treat illegal character
c
150 buff(i) = '0'
valr8 = 0.0d0
icond = 1
c
1000 if( dpoint )then
valr8 = valr8/(ten**ip)
endif
c
if( .not.plus )then
valr8 = -valr8
endif
c
return
end
subroutine todmsp(val,id,im,s,isign)
*** convert position radians to deg,min,sec
*** range is [-pi to +pi]
implicit double precision(a-h,o-z)
common/const/pi,rad
1 if(val.gt.pi) then
val=val-pi-pi
go to 1
endif
2 if(val.lt.-pi) then
val=val+pi+pi
go to 2
endif
if(val.lt.0.d0) then
isign=-1
else
isign=+1
endif
s=dabs(val*rad)
id=idint(s)
s=(s-id)*60.d0
im=idint(s)
s=(s-im)*60.d0
*** account for rounding error
is=idnint(s*1.d5)
if(is.ge.6000000) then
s=0.d0
im=im+1
endif
if(im.ge.60) then
im=0
id=id+1
endif
return
end
subroutine trim (buff,lgh,hem)
c
implicit integer (i-n)
character*1 ch,hem
character*1 buff(*)
integer*4 lgh
c
ibeg = 1
do 10 i=1,50
if( buff(i).ne.' ' )then
goto 11
endif
ibeg = ibeg+1
10 continue
11 continue
if( ibeg.ge.50 )then
ibeg = 1
buff(ibeg) = '0'
endif
c
iend = 50
do 20 i=1,50
j = 51-i
if( buff(j).eq.' ' )then
iend = iend-1
else
goto 21
endif
20 continue
21 continue
c
ch = buff(ibeg)
if( hem.eq.'N' )then
if( ch.eq.'N' .or. ch.eq.'n' .or. ch.eq.'+' )then
hem = 'N'
ibeg = ibeg+1
endif
if( ch.eq.'S' .or. ch.eq.'s' .or. ch.eq.'-' )then
hem = 'S'
ibeg = ibeg+1
endif
c
c check for wrong hemisphere entry
c
if( ch.eq.'E' .or. ch.eq.'e' )then
hem = '*'
ibeg = ibeg+1
endif
if( ch.eq.'W' .or. ch.eq.'w' )then
hem = '*'
ibeg = ibeg+1
endif
elseif( hem.eq.'W' )then
if( ch.eq.'E' .or. ch.eq.'e' .or. ch.eq.'+' )then
hem = 'E'
ibeg = ibeg+1
endif
if( ch.eq.'W' .or. ch.eq.'w' .or. ch.eq.'-' )then
hem = 'W'
ibeg = ibeg+1
endif
c
c check for wrong hemisphere entry
c
if( ch.eq.'N' .or. ch.eq.'n' )then
hem = '*'
ibeg = ibeg+1
endif
if( ch.eq.'S' .or. ch.eq.'s' )then
hem = '*'
ibeg = ibeg+1
endif
elseif( hem.eq.'A' )then
if( .not.('0'.le.ch .and. ch.le.'9') )then
hem = '*'
ibeg = ibeg+1
endif
else
c do nothing
endif
c
c
do 30 i=ibeg,iend
ch = buff(i)
c
if( ch.eq.':' .or. ch.eq.'.' )then
goto 30
elseif( ch.eq.' ' .or. ch.eq.',' )then
buff(i) = ':'
elseif( '0'.le.ch .and. ch.le.'9' )then
goto 30
else
buff(i) = ':'
endif
c
30 continue
c
c left justify buff() array to its first character position
c also check for a ":" char in the starting position,
c if found!! skip it
c
j = 0
ib = ibeg
ie = iend
c
do 40 i=ib,ie
if( i.eq.ibeg .and. buff(i).eq.':' )then
c
c move the 1st position pointer to the next char &
c do not put ":" char in buff(j) array where j=1
c
ibeg = ibeg+1
goto 40
endif
j = j+1
buff(j) = buff(i)
40 continue
c
c
lgh = iend-ibeg+1
j = lgh+1
buff(j) = '$'
c
c clean-up the rest of the buff() array
c
do 50 i=j+1,50
buff(i) = ' '
50 continue
c
c save a maximum of 20 characters
c
if( lgh.gt.20 )then
lgh = 20
j = lgh+1
buff(j) = '$'
endif
c
return
end