SUBROUTINE STGERP(A,B,PHI,PHIO,CHI,CHIO,DSLON,CM,DM,ICODE,KO,C1,C2        
     1                  ,R)                                                     
C                                                                               
C                                                                               
C        THIS ROUTINE COMPUTES THE TWO BY TWO COVARIANCE MATRIX OF THE          
C      STEREOGRAPHIC (DOUBLE PROJECTION) PLANE COORDINATES GIVEN THE            
C      CORRESPONDING TWO BY TWO COVARIANCE MATRIX OF THE ELLIPSOIDAL            
C      COORDINATES OF THE POINT; AND CONVERSELY.                                
C                                                                               
C        INPUT:                                                                 
C                                                                               
C              A,B   - SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE          
C                      ELLIPSOID, IN METRES.                                    
C              PHI   - ELLIPSOIDAL LATITUDE OF THE POINT, IN RADIANS.           
C              PHIO  - ELLIPSOIDAL LATITUDE OF THE ORIGIN OF THE                
C                      PROJECTION, IN RADIANS.                                  
C              CHI   - SPHERICAL LATITUDE OF THE POINT, IN RADIANS.             
C              CHIO  - SPHERICAL LATITUDE OF THE ORIGIN OF THE PROJECTION,      
C                      IN RADIANS.                                              
C              DSLON - SPHERICAL LONGITUDE OF THE POINT MINUS THE               
C                      SPHERICAL LONGITUDE OF THE ORIGIN OF THE                 
C                      PROJECTION, IN RADIANS. (LONGITUDE IS MEASURED           
C                      POSITIVE EAST).                                          
C              CM    - COVARIANCE MATRIX (2X2) OF THE PLANE STEREOGRAPHIC       
C                      GRID COORDINATES, IN METRES SQUARED. (SEE ICODE          
C                      BELOW).                                                  
C              DM    - COVARIANCE MATRIX (2X2) OF THE ELLIPSOIDAL               
C                      COORDINATES, IN RADIANSSQUARED. (SEE ICODE).             
C              ICODE - MAY BE EITHER 1 OR -1.                                   
C                    - WHEN ICODE IS 1 CM WILL BE COMPUTED GIVEN DM.            
C                      (DM IS INPUT AND CM IS OUTPUT).                          
C                    - WHEN ICODE IS -1 DM WILL BE COMPUTED GIVEN CM.           
C              KO    - POINT SCALE FACTOR AT THE ORIGIN.                        
C              C1    - A CONSTANT COMPUTED IN STGINL.                           
C              C2    - A CONSTANT COMPUTED IN STGINL.                           
C              R     - RADIUS OF THE CONFORMAL SPHERE,(COMPUTED IN STGINL).     
C                                                                               
C        OUTPUT:                                                                
C                                                                               
C              CM OR DM DEPENDING ON THE VALUE OF ICODE; (1 OR -1 )             
C                                                                               
C        THIS ROUTINE UTILIZES RIGOROUS EQUATIONS FOR THE ERROR                 
C      PROPAGATION; APPROXIMATIONS MAY BE MORE EFFICIENT AND YEILD              
C      SUFFICIENT ACCURACY IN THE RESULTS.                                      
C                                                                               
C                                      WRITTEN BY R.R.STEEVES,JULY,1977.        
C                                                                               
      IMPLICIT REAL* 8 (A-H,J-Z)                                                
      DIMENSION CM(2,2),DM(2,2)                                                 
      E2=(A*A-B*B)/(A*A)                                                        
      E=DSQRT(E2)                                                               
      ESP=E*DSIN(PHI)                                                           
      RR=3.141592653589793D0/4.D0+PHI/2.D0                                      
      Q=((1.D0-ESP)/(1.D0+ESP))**(E/2.D0)                                       
      R2=2.D0*KO*R                                                              
      SC=DSIN(CHI)                                                              
      SCO=DSIN(CHIO)                                                            
      CC=DCOS(CHI)                                                              
      CCO=DCOS(CHIO)                                                            
      CDS=DCOS(DSLON)                                                           
      M=R2*DSIN (DSLON)*(SC+SCO)                                                
      N=(1.D0+SC*SCO+CC*CCO*CDS)**2                                             
      P=R2*(CDS*(1.D0+SC*SCO)+CC*CCO)                                           
      TR=DTAN(RR)                                                               
      U=C1*C2*(Q*TR)**(C1-1.D0)*Q*(1.D0/2.D0/DCOS(RR)**2-E2*DCOS(PHI)*TR        
     1/(1.D0-E2*DSIN(PHI)**2))/(1.D0+C2**2*(Q*TR)**(2.D0*C1))                   
      U=U*2.D0                                                                  
      CNC=C1/N*CC                                                               
      J11=-M/N*U                                                                
      J12=CNC*P                                                                 
      J21=P/N*U                                                                 
      J22=CNC*M                                                                 
      IF(ICODE.LT.0)GOTO1                                                       
      CM(1,1)=J11**2*DM(1,1)+2.D0*J11*J12*DM(1,2)+J12**2*DM(2,2)                
      CM(1,2)=J11*J21*DM(1,1)+J12*J21*DM(1,2)+J11*J22*DM(1,2)+J12*J22*DM        
     1        (2,2)                                                             
      CM(2,1)=CM(1,2)                                                           
      CM(2,2)=J21**2*DM(1,1)+2.D0*J21*J22*DM(1,2)+J22**2*DM(2,2)                
      GO TO 10                                                                  
    1 CONTINUE                                                                  
      DET=J11*J22-J21*J12                                                       
      J1=J11                                                                    
      J11=J22/DET                                                               
      J12=-J12/DET                                                              
      J21=-J21/DET                                                              
      J22=J1/DET                                                                
      DM(1,1)=J11**2*CM(1,1)+2.D0*J11*J12*CM(1,2)+J12**2*CM(2,2)                
      DM(1,2)=J11*J21*CM(1,1)+J12*J21*CM(1,2)+J11*J22*CM(1,2)+J12*J22*CM        
     1        (2,2)                                                             
      DM(2,1)=DM(1,2)                                                           
      DM(2,2)=J21**2*CM(1,1)+2.D0*J21*J22*CM(1,2)+J22**2*CM(2,2)                
   10 RETURN                                                                    
      END