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