SUBROUTINE ELLINV(PHI1,AMB1,PHI2,AMB2,ALP,DIS,ALP2,A,B)                   
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      REAL*8 M1,N1,N2                                                           
C                                                                               
C     SUBROUTINE  TO DO INVERSE PROBLEM ON THE ELLIPSOID USING                  
C      PUISSANTS FORMULAE                                                       
C                                                                               
C     ***INPUT***                                                               
C                                                                               
C     PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS)                               
C     AMB1----GEODETIC LONGITUDE OF PT. 1 (RADIANS)                             
C     PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS)                               
C     AMB2----GEODETIC LONGITUDE OF PT. 2 (RADIANS)                             
C     A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES)                       
C     B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES)                       
C                                                                               
C     ***OUTPUT***                                                              
C                                                                               
C     DIS----ELLIPSOID DISTANCE(METRES)                                         
C     ALP----GEODETIC AZIMUTH FROM PT.1 TO PT.2 (RADIANS)                       
C     ALP2----COMPUTED GEODETIC AZIMUTH FROM PT.2 TO PT.1 (RADIANS)             
C                                                                               
      WRITE(6,5)                                                                
  5    FORMAT('0',20X,'ITERATIONS IN PUISSANTS INVERSE PROBLEM')                
      PI=DARCOS(-1.0D0)                                                         
      RO=6.48D5/PI                                                              
      E=(A**2-B**2)/(A**2)                                                      
      N1=A/((1.0D0-E*((DSIN(PHI1))**2))**.5D0)                                  
      M1=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI1))**2))**1.5D0)                       
      N2=A/((1.0D0-E*((DSIN(PHI2))**2))**.5D0)                                  
      DP=PHI2-PHI1                                                              
      DL=AMB2-AMB1                                                              
      F1=DP*M1/(1.0D0-3.0D0*E*DSIN(PHI1)*DCOS(PHI1)*DP/(2.0D0*(1.0D0            
     D-E*((DSIN(PHI1))**2))))                                                   
      F2=DL*N2*DCOS(PHI2)                                                       
      ALP=DATAN2(F2,F1)                                                         
      IF(ALP.LT.0.0D0)ALP=2.0D0*PI+ALP                                          
      AAA=ALP                                                                   
      DIS=F1/DCOS(ALP)                                                          
      DO 1 JJ=1,10                                                              
      F11=(DIS**2)*DTAN(PHI1)*((DSIN(ALP))**2)/(2.0D0*N1)                       
      F111=(DIS**3)*DCOS(ALP)*((DSIN(ALP))**2)*(1.0D0+3.0D0*((DTAN(             
     NPHI1))**2))/(6.0D0*(N1**2))                                               
      F22=(DIS**3)*DSIN(ALP)/(6.0D0*(N2**2))                                    
      F222=(DIS**3)*((DSIN(ALP))**3)/(6.0D0*(N2**2)*((DCOS(PHI2))**2))          
      F=F1+F11+F111                                                             
      FP=F2+F22-F222                                                            
      ALP=DATAN2(FP,F)                                                          
      IF(ALP.LT.0.0D0)ALP=2.0D0*PI+ALP                                          
      DIS=F/DCOS(ALP)                                                           
      DIFF=ALP-AAA                                                              
       WRITE(6,3)DIFF                                                           
  3    FORMAT('0',20X,'AZIMUTH DIFFERENCE=',F12.9)                              
       CALL RADMS(DIFF,ID,IM,SE)                                                
        WRITE(6,4)ID,IM,SE                                                      
  4   FORMAT('0',20X,'AZIMUTH DIFFERENCE=',2X,I4,2X,I3,2X,F9.5)                 
       CALL RADMS(ALP,IDEG,IMIN,SEC)                                            
       WRITE(6,6)IDEG,IMIN,SEC                                                  
  6     FORMAT('0',20X,'AZIMUTH =',2X,I4,2X,I3,2X,F7.3)                         
      WRITE(6,7)DIS                                                             
  7    FORMAT('0',20X,'ELLIPSOID DISTANCE=',F15.4)                              
      IF(DABS(DIFF).LT.1.0D-9)GO TO 10                                          
      AAA=ALP                                                                   
   1  CONTINUE                                                                  
  10  CONTINUE                                                                  
      PM=(PHI1+PHI2)/2.0D0                                                      
      DP=DP*RO                                                                  
      DL=DL*RO                                                                  
      DA1=DL/RO*DSIN(PM)/(DCOS(DP/(2.0D0*RO)))                                  
      DA2=((DL/RO)**3)/12.0D0*(DSIN(PM)/(DCOS(DP/(2.0D0*RO)))                   
     N-((DSIN(PM))**3)/((DCOS(DP/(2.0D0*RO)))**3))                              
      DA=(DA1+DA2)*RO                                                           
      ALP2=ALP+PI+DA/RO                                                         
      IF(ALP.GT.PI)ALP2=ALP-PI+DA/RO                                            
      IF(JJ.GT.9)WRITE(6,100)                                                   
  100 FORMAT('0','WARNING---MAXIMUM # OF ITERATIONS PERFORMED IN ELLINV         
     CSUBROUTINE---CHECK FOR ERRORS')                                           
      RETURN                                                                    
      END