PRO ipolster,alon,alat,x,y,xlam1,slat1 ; ; COMPUTES THE INVERSE POLAR STEROGRAPHIC TRANSFORMATION FOR (X,Y) ; GIVEN IN KM WITH REFERENCES LON,LAT=(XLAM,SLAT). ; OUTPUT LON,LAT=ALON,ALAT ; ; ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES ; AS RECEIVED FROM M. DRINKWATER (1994). UPDATED BY D. LONG TO ; IMPROVE ACCURACY USING ITERATION WITH FORWARD TRANSFORM. ; ; version 1.1 modified - M. Drinkwater - JPL 9/23/97 ; corrected residual fortran "SIGN" function calls ; version 1.2 modified - D. Long - BYU 4/30/02 ; added DOUBLE() to inputs to ensure correct values in new IDL versions ; E2=0.006693883D0 RE=6378.273D0 PI2=1.570796327D0 DTR=PI2/90.0D0 ; ; FIRST USE APPROXIMATE INVERSE CALCULATION ; E=SQRT(E2) E22=E2*E2 E23=E2*E2*E2 X1=DOUBLE(X) Y1=DOUBLE(Y) XLAM=DOUBLE(XLAM1) SLAT=DOUBLE(SLAT1) RHO=X1*X1+Y1*Y1 IF RHO GT 0.0 THEN RHO=SQRT(RHO) IF RHO LT 0.05 THEN BEGIN ALON=XLAM ; ALAT=SIGN(90.0,SLAT) ; this is a fortran statement IF(SLAT GE 0.) THEN ALAT=90.0 IF(SLAT LT 0.) THEN ALAT=-90.0 ENDIF ELSE BEGIN SN=1.0D0 SLAT1=SLAT IF SLAT LT 0.0 THEN BEGIN SN=-1.0D0 SLAT1=-SLAT ENDIF CM=COS(SLAT1 * DTR)/SQRT(1.0D0-E2*SIN(SLAT1 * DTR)^2) T=TAN(DTR*(45.0D0-0.5*SLAT1))/ $ ((1.D0-E*SIN(SLAT1*DTR))/(1.D0+E*SIN(SLAT1*DTR)))^(E*0.5D0) T=RHO*T/(RE*CM) CHI=PI2-2.0D0*ATAN(T) T=CHI+ $ (0.5D0*E2+5.0D0*E22/24.0D0+E23/12.0D0)*SIN(2.0D0*CHI)+ $ (7.0D0*E22/48.0D0+29.0D0*E23/240.0D0)*SIN(4.0D0*CHI)+ $ (7.0D0*E23/120.0D0)*SIN(6.0D0*CHI) ALAT=SN*(T*90.0D0/PI2) ALON=SN*ATAN(SN*X1,-SN*Y1)/DTR+XLAM IF ALON LT -180.0 THEN ALON=ALON+360.0 IF ALON GT 180.0 THEN ALON=ALON-360.0 ENDELSE ; ; USING THE APPROXIMATE RESULT AS A STARTING POINT, ITERATE TO IMPROVE ; THE ACCURACY OF THE INVERSE SOLUTION ; SN1=1.0 IF SLAT LT 0.0 THEN SN1=-1.0 A=ATAN(Y1,X1) / DTR R=SQRT(X1*X1+Y1*Y1) ICNT=0 S: polster,ALON,ALAT,XX,YY,XLAM,SLAT RR=SQRT(XX*XX+YY*YY) RERR=SN1*(RR-R)/180.0D0 AA=ATAN(YY,XX)/DTR AERR=AA-A IF ABS(AERR) GT 180.0 THEN AERR=360.0-AERR ; ; CHECK FOR CONVERGENCE ; IF (((ABS(RERR) LT 0.001) AND (ABS(AERR) LT 0.001)) OR (ICNT GT 9)) THEN GOTO, S40 ; ; CONSTRAIN UPDATES ; ALON=ALON+AERR IF ABS(ALON) GT 360.0 THEN ALON= ALON MOD 360.0 IF ALAT*SLAT LT 0.0 THEN BEGIN RERR=RERR*(1.0D0-SIN(DTR*ABS(ALAT))) ; IF ABS(RERR) GT 2.0 THEN RERR=SIGN(2.0,RERR)/ICNT IF ABS(RERR) GT 2.0 THEN BEGIN IF(RERR GE 0.) THEN RERR=2.0/FLOAT(ICNT) IF(RERR LT 0.) THEN RERR=-2.0/FLOAT(ICNT) ENDIF ENDIF ALAT=ALAT+RERR IF ABS(ALAT) GT 90.0 THEN ALAT=SIGN(90.0,ALAT) ICNT=ICNT+1 GOTO, S S40: IF ABS(ALON) GT 360.0 THEN ALON=ALON MOD 360.0 END