cb::forward c program forward c c********1*********2*********3*********4*********5*********6*********7** c c name: forward c version: 200208.19 c author: stephen j. frakes c purpose: to compute a geodetic forward (direct problem) 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 arc meridional arc distance latitude p1 to p2 (meters) c b semiminor axis polar (in meters) c baz azimuth back (in radians) c blimit geodetic distance allowed on ellipsoid (in meters) c buff input char buffer array c dd,dm,ds temporary values for degrees, minutes, seconds c dmt,d_mt char constants for units (in meters) c dd_max maximum ellipsoid distance -1 (in meters) 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 value for ellipsoid distance 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 azd1,azm1,saz1 forward azimuth - degrees,minutes,seconds c iazsn forward azimuth - flag (+/- 1) c glat1,glon1 station one - (lat & lon in radians ) c c name2 name of station two c ld2,lm2,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 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 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, dirct1, todmsp c gethem, trim, bufdms, gvalr8, gvali4, fixdms, gpnarc 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:: at request of Dave Doyle c::200207.18, rws, modified i/o & standardized program documentation c:: added subs trim, bufdms, gethem, gvali4, gvalr8 c::200207.23, rws, added sub gpnarc 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 c:: - added logical nowebb c::200208.xx, sjf, program version number 2.0 c********1*********2*********3*********4*********5*********6*********7** ce::forward 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 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 5 write(*,*) ' Program Forward - 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 r1 = 0.0d0 r2 = pi/2.0d0 call gpnarc ( a, f, esq, pi, r1, r2, arc ) c c compute the geodetic limit distance (blimit), it is equal c to twice the distance to the pole minus one meter c blimit = 2.0d0*arc-1.0d0 c c maximum distance allowed on ellipsoid c dd_max = blimit 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 )' c 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 20 write(*,*) 'DDD MM SS.sss Forward Azimuth : (from north)' c 24 hem='A' read(*,11) buff call trim (buff,lgh,hem) call bufdms (buff,lgh,hem,dd,dm,ds,ierror) c if( ierror.eq.0 )then iazsn = 1 irazi1 = 0 else irazi1 = 1 write(*,*) ' Invalid Azimuth ... Please re-enter it ' write(*,*) ' ' if( nowebb )then goto 24 endif endif c azd1 = dd azm1 = dm saz1 = ds c call getrad(azd1,azm1,saz1,iazsn,faz) c write(*,*) 'DDDDDD.dddd Ellipsoidal Distance : (in meters)' c 25 ss = 0.0d0 read(*,*) ss ss = dabs(ss) c if( ss.lt.dd_max )then edist = ss irdst1 = 0 else irdst1 = 1 write(*,*) ' Invalid Distance ... Please re-enter it ' write(*,*) ' ' if( nowebb )then goto 25 endif edist = 0.001d0 endif c call dirct1 (glat1,glon1,glat2,glon2,faz,baz,edist) c c set the azimuth seconds tolerance 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 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 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: Forward Azimuth .. Invalid Entry ') 87 format(' Error: Ellipsoid Distance .. Invalid Entry ') 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...) ') 93 format(' Geodetic distance is too long ') c write(3,*) ' ' write(3,41) name1 write(3,*) ' ---------------- ' c 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,*) ' ---------------- ' write(3,18) ld2, lm2, sl2, d_ns2 write(3,19) lod2,lom2,slo2,d_ew2 c write(3,*) ' ' if( irazi1.eq.0 )then write(3,34) iaz1,maz1,saz1 else write(3,88) write(3,86) write(3,92) write(3,88) endif c write(3,35) iaz2,maz2,saz2 c if( irdst1.eq.0 )then write(3,32) edist,dmt else write(3,88) write(3,87) write(3,93) write(3,88) endif c write(3,*) ' ' endif c write(*,*) ' ' write(*,*) ' ' write(*,*) ' ' write(*,*) ' 1) Another forward, different ellipsoid.' write(*,*) ' 2) Same ellipsoid, two new stations.' write(*,*) ' 3) Same ellipsoid, same First Station.' write(*,*) ' 4) Let the Second Station be the First Station.' write(*,*) ' 5) 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 elseif(answer.eq.'4')then glat1 = glat2 glon1 = glon2 goto 20 else write(*,*) ' Thats all, folks!' endif c 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 DIRCT1(GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S) C C *** SOLUTION OF THE GEODETIC DIRECT 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 ANTIPODAL C C *** A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID C *** F IS THE FLATTENING OF THE REFERENCE ELLIPSOID C *** LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST C *** AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH C *** GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A C C *** PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75 C *** MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608 C IMPLICIT REAL*8 (A-H,O-Z) COMMON/CONST/PI,RAD COMMON/ELIPSOID/A,F DATA EPS/0.5D-13/ R=1.-F TU=R*DSIN(GLAT1)/DCOS(GLAT1) SF=DSIN(FAZ) CF=DCOS(FAZ) BAZ=0. IF(CF.NE.0.) BAZ=DATAN2(TU,CF)*2. CU=1./DSQRT(TU*TU+1.) SU=TU*CU SA=CU*SF C2A=-SA*SA+1. X=DSQRT((1./R/R-1.)*C2A+1.)+1. X=(X-2.)/X C=1.-X C=(X*X/4.+1)/C D=(0.375D0*X*X-1.)*X TU=S/R/A/C Y=TU 100 SY=DSIN(Y) CY=DCOS(Y) CZ=DCOS(BAZ+Y) E=CZ*CZ*2.-1. C=Y X=E*CY Y=E+E-1. Y=(((SY*SY*4.-3.)*Y*CZ*D/6.+X)*D/4.-CZ)*SY*D+TU IF(DABS(Y-C).GT.EPS)GO TO 100 BAZ=CU*CY*CF-SU*SY C=R*DSQRT(SA*SA+BAZ*BAZ) D=SU*CY+CU*SY*CF GLAT2=DATAN2(D,C) C=CU*CY-SU*SY*CF X=DATAN2(SY*SF,C) C=((-3.*C2A+4.)*F+4.)*C2A*F/16. D=((E*CY*C+CZ)*SY*C+Y)*SA GLON2=GLON1+X-(1.-C)*D*F BAZ=DATAN2(SA,BAZ)+PI 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 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) 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 real*8 function zvalue ( ss, tol ) c implicit double precision (a-z) c c to check for dx,dy,dz or c dn,de,du values below tol. c if( dabs(ss).lt.tol )then ss = 0.0d0 endif c zvalue = ss c return end