C SUBROUTINE OMXYPL(X,Y,A,B,C,KC,AELLIP,BELLIP,ALPHA0,LAMDA0, &ALPHAC,V,U,PHI,ELAMDA) C C THIS SUBROUTINE COMPUTES THE SKEW COORDINATES V,U AND THE C GEODETIC COORDINATES PHI AND LAMDA GIVEN THE CARTESIAN C COORDINATES X,Y OF THE OBLIQUE MERCATOR PROJECTION. C THE FORMULAE ARE BASED ON THE HOTINE'S METHOD OF APOSPHERIC C PROJECTION (1946-47). C C INPUT: C X,Y--CARTESIAN COORDINATES C A,B,C--APOSPHERIC CONSTANTS C KC--SCALE FACTOR AT THE CENTRE OF PROJECTION C ALPHA0--ANGLE,AT THE ORIGIN,BETWEEN THE MERIDIAN LAMDA0 C AND THE AXIS OF PROJECTION C AELLIP,BELLIP--SEMI-MAJOR AND SEMI-MINOR AXIS OF THE C ELLIPSOID RESPECTIVELY C LAMDA0--MERIDIAN THAT PASSES THROUGH THE ORIGIN OF THE C PROJECTION C ALPHAC--ANGLE(IN RADIANS),AT THE PROJECTION CENTRE,WHICH C THE MERIDIAN MAKES WITH THE AXIS OF PROJECTION. C OUTPUT: C V,U--SKEW GRID COORDINATES C PHI--GEODETIC LATITUDE OF THE POINT,IN RADIANS C ELAMDA--GEODETIC LONGITUDE(POSITIVE EAST OF GREENWICH) C OF THE POINT,IN RADIANS C C WRITTEN BY A.L.KOK,AUGUST,1982 C IMPLICIT REAL*8(A-H,K-Z) C PI=DARCOS(-1.0D0) F=DSIN(ALPHA0) G=DCOS(ALPHA0) C TAC=DTAN(ALPHAC) SAC=DSIN(ALPHAC) CAC=DCOS(ALPHAC) C D=A*KC/B E=DSQRT((AELLIP**2-BELLIP**2)/AELLIP**2) C U= X*SAC+Y*CAC V= X*CAC-Y*SAC C P=DSIN(U/D) R=DSINH(V/D) S=DCOSH(V/D) C TEMP1=S+G*P-F*R TEMP2=S-G*P+F*R TEMP3=DLOG(TEMP1/TEMP2) PSI=(1.0D0/(2.0D0*B))*TEMP3-(C/B) C C NEWTON-RAPHSON ITERATION METHOD USED TO COMPUTE LATITUDE C EQ=DEXP(PSI) PHI=2.0D0*DATAN(EQ)-PI/2.0D0 100 SF=DSIN(PHI) F1=0.5D0*(DLOG(1+SF)-DLOG(1-SF)+E*DLOG(1-E*SF)-E*DLOG(1+E*SF))-PSI E2=E**2 F2=(1-E2)/((1-E2*SF**2)*DCOS(PHI)) DELPHI=F1/F2 PHI=PHI-DELPHI IF (DABS(DELPHI) .GT. 1.0D-12) GO TO 100 C C TO COMPUTE THE LONGITUDE TEMP1= -G*R-F*P TEMP2=(-1.0D0)*TEMP1/DCOS(U/D) ELAMDA=LAMDA0+(1.0D0/B)*DATAN(TEMP2) C RETURN END C