PRO GEOTRANS,BASIN,LAT,LON,ICODE,X,Y

;  From norton@glerl.noaa.gov Thu Jan 30 10:33:51 EST 1997
;  fortran converted to IDL, 8/2001, G. Lang

;      SUBROUTINE FOR GEOGRAPHIC TRANSFORMATION   
; 
;      AN ERIM SUBROUTINE MODIFIED FOR GLERL BY DAVID NORTON              
; 
;      THIS PROGRAM CAN BE USED TO FIND THE DISTANCE FROM ONE POINT TO    
;      ANOTHER IN THE GREAT LAKES REGION, OR THE LATITUDE AND LONGITUDE   
;      OF A POINT IF ITS DISTANCE FROM A KNOWN LATITUDE AND LONGITUDE     
;      IS KNOWN.  THE ORIGIN MUST BE TO THE NORTH WEST OF THE OTHER       
;      POINTS.  IT REQUIRES DECIMAL DEGREES AND DISTANCES IN METERS.      
;      THE FIRST CALL MUST BE FOR THE ORIGIN.  EXAMPLE:                   
; 
;      BASIN = 1.     (2. FOR MICH, 3. FOR HURON, ETC.)                   
;      IF A MAP BASE OTHER THAN DAVID NORTON'S IS USED, MAKE BASIN        
;      EQUAL TO THE "CENTRAL MERIDIAN" OF THE MAP (EX. -83.5).            
;      ICODE = 1          
;      FLT =              
;      FLN =      (NEGATIVE AT THIS LONGITUDE)    
;      CALL GT(BASIN,FLT,FLN,ICODE,XFO,YFO)       
; 
;      SINCE XFO AND YFO ARE DISTANCES FROM THE EQUATOR, THE SUBSEQUENT   
;      CALLS FOR LT,LN OR X,Y MUST BE ADJUSTED AS FOLLOWS:                
; 
;      IF YOU DESIRE A DISTANCE FROM THE ORIGIN, MAIN PROGRAM CODING      
;      SHOULD LOOK LIKE:  
; 
;      REAL LT,LN         
;      LT =               
;      LN =        (NEGATIVE AT THIS LONGITUDE)   
;      CALL GT(BASIN,LT,LN,ICODE,X,Y)             
;      X = X-XFO          
;      Y = YFO-Y          
; 
;      IF YOU DESIRE A LATITUDE, LONGITUDE, MAIN PROGRAM CODING           
;      SHOULD LOOK LIKE:  
; 
;      REAL LT,LN         
;      ICODE = (-1)       
;      X =                
;      Y =                
;      X = XFO+X          
;      Y = YFO-Y          
;      CALL GT(BASIN,LT,LN,ICODE,X,Y)             

PCM=FLTARR(6)
OLAT=FLTARR(6)
OLON=FLTARR(6)

RS=6371100.
DPR=57.29578
PO4=.785398
A15=6378206.4
B1=6356583.8
PCM=[-88.0,-87.0,-82.0,-82.0,-82.0,-77.5]

;     CENTRAL MERIDIANS PCM/SUPERIOR,MICH.,HURON,ST.CLAIR,ERIE,ONT./     

OLAT=[51.0,47.0,48.0,44.0,44.5,45.6667]
OLON=[-93.75,-90.25,-85.6667,-84.0,-85.75,-80.5]

;     ORIGIN LATITUDES AND LONGITUDES FOR EACH LAKE.                     
;     ICODE IS + TO GET X,Y OR - FOR LON,LOT       
;
      
IF BASIN EQ 0 THEN begin                 
   PCM(0) = BASIN                     
   LAKE = 1           
endif else begin
   LAKE = BASIN         
   IF LON eq 0.0 THEN begin
      LAT=OLAT(LAKE-1)      
      LON=OLON(LAKE-1)      
   end
endelse         

;     CHECK FOR X,Y OR LON,LAT DESIRED  

IF ICODE GT 0 then begin                  

;     COMPUTE X,Y        

   P=LAT/DPR           
   TP=TAN(P)          
   SNPL=SIN(P)*(LON-PCM(LAKE-1))/DPR             
   X=RS*SIN(SNPL)/TP  
   Y=RS*(P+(1-COS(SNPL))/TP)                  
   RETURN             

end       

;     COMPUTE LAT,LON      

LAT=Y/RS            
LON=X/(RS*COS(LAT))  

for I=1,10 do begin                      

   TP=TAN(LAT)         
   SNPL=SIN(LAT)*LON    
   XN=RS*SIN(SNPL)/TP 
   YN=RS*(LAT+(1.-COS(SNPL))/TP)               
   DY=YN-Y            
   DX=XN-X            
   LAT=LAT-DY/RS        
   LON=LON-DX/(COS(LAT)*RS)                      

;     CHECK FOR CONVERGENCE                      

   IF DX LT 1.0 AND DY LT 1.0 then begin
      LAT=LAT*DPR          
      LON=LON*DPR+PCM(LAKE-1)
      RETURN             
   end
end          

print, '  WARNING...  CONVERGENCE ERRORS OF ',DX,DY 

return
end      

; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; Main Program
; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

; Reads pairs of lat,lon and prints dist from previous pair and cum dist (km)

ICODE = 1          
BASIN = 1.

filein='mapdist.in'
spawn,'wc -l '+filein,a
reads,a,nrec
a=fltarr(2,nrec)
openr,1,filein
readf,1,a
close,1
a=transpose(a)
lat=a(*,0)
lon=a(*,1)

dist=0.
cumdist=0.
delX=0.
delY=0.
for i=0,nrec-1 do begin
   GEOTRANS,BASIN,LAT(i),LON(i),ICODE,X,Y
   if i gt 0 then begin
      delx = x-xprev         
      dely = yprev-y          
      dist=sqrt(delx*delx+dely*dely)
   end
   xprev=x
   yprev=y
   cumdist=cumdist+dist
   print,LAT(i),LON(i),delX/1000.,delY/1000.,dist/1000.,cumdist/1000., $
   format="(2f10.4,2f8.2,',  Dist (km) =',f8.2,',  Cum Dist (km) =',f8.2)"
end

end