SUBROUTINE ELLDIR(PHI,AMB,ALP,DIS,ALP2,PHI2,AMB2,A,B)                     
C                                                                               
C     SUBROUTINE TO DO DIRECT PROBLEM ON THE ELLIPSOID USING                    
C      PUISSANTS FORMULAE                                                       
C                                                                               
C     ***INPUT***                                                               
C     PHI----GEODETIC LATITUDE OF PT.1 (RADIANS)                                
C     AMB----GEODETIC LONGITUDE OF PT. 1 (RADIANS)                              
C     DIS----GIVEN ELLIPSOID DISTANCE(METRES)                                   
C     ALP----GIVEN GEODETIC AZIMUTH FROM PT.1 TO 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     PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS)                               
C     AMB2----GEODETIC LONGITUDE OF PT. 2 (RADIANS)                             
C     ALP2----COMPUTED GEODETIC AZIMUTH FROM PT.2 TO PT.1 (RADIANS)             
C                                                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      REAL*8 N1,N2,M1                                                           
       WRITE(6,5)                                                               
      E=(A**2-B**2)/(A**2)                                                      
      PI=DARCOS(-1.0D0)                                                         
      RO=6.48D5/PI                                                              
      N1=A/((1.0D0-E*((DSIN(PHI))**2))**.5D0)                                   
      M1=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI))**2))**1.5D0)                        
      DP1=DIS/N1*DCOS(ALP)                                                      
      DP2=(DIS**2)*DTAN(PHI)*((DSIN(ALP))**2)/(2.0D0*(N1**2))                   
      DP3=(DIS**3)*DCOS(ALP)*((DSIN(ALP))**2)*(1.0D0+((DTAN(PHI))**2)           
     C*3.0D0)/(6.0D0*(N1**3.0D0))                                               
      DP=(DP1-DP2-DP3)*RO                                                       
      DO 1 JJ=1,10                                                              
      DP1=DIS*DCOS(ALP)/M1                                                      
      DP2=(DIS**2)*DTAN(PHI)*((DSIN(ALP))**2)/(2.0D0*M1*N1)                     
      DP3=(DIS**3)*DCOS(ALP)*((DSIN(ALP))**2)*(1.0D0+3.0D0*                     
     F((DTAN(PHI))**2))/(6.0D0*M1*(N1**2))                                      
      DS=RO*(DP1-DP2-DP3)                                                       
      DP4=(1.0D0-3.0D0*E*DSIN(PHI)*DCOS(PHI)/(2.0D0*(1.0D0-E*                   
     ^((DSIN(PHI))**2)))*DP/RO)                                                 
      DS=DS*DP4                                                                 
      PHI2=PHI+DS/RO                                                            
      DIFF=DP-DS                                                                
        DIFF =DIFF/RO                                                           
       WRITE(6,3)DIFF                                                           
       CALL RADMS(DIFF,ID,IM,SE)                                                
      WRITE(6,4)ID,IM,SE                                                        
  4     FORMAT('0',20X,'DELTA PHI =',2X,I3,2X,I3,2X,F9.5)                       
        DIFF=DIFF*RO                                                            
  3    FORMAT('0',20X,'DELTA PHI=',F12.9)                                       
  5     FORMAT('0',20X,'ITERATIONS IN PUISSANTS DIRECT PROBLEM')                
      IF(DABS(DIFF).LT.1.0D-5)GO TO 2                                           
      DP=DS                                                                     
   1  CONTINUE                                                                  
  2   CONTINUE                                                                  
      N2=A/((1.0D0-E*((DSIN(PHI2))**2))**.5D0)                                  
      DL1=DIS*DSIN(ALP)/(N2*DCOS(PHI2))                                         
      DL2=1.0D0-(DIS**2)*(1.0D0-((DSIN(ALP))**2)/((DCOS(PHI2))**2))             
     B/(6.0D0*(N2**2))                                                          
      DL=RO*DL1*DL2                                                             
      AMB2=AMB+DL/RO                                                            
      PM=(PHI+PHI2)/2.0D0                                                       
      DA1=DL/RO*DSIN(PM)/(DCOS(DS/(2.0D0*RO)))                                  
      DA2=((DL/RO)**3)/12.0D0*(DSIN(PM)/(DCOS(DS/(2.0D0*RO)))                   
     N-((DSIN(PM))**3)/((DCOS(DS/(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 ELLDIR         
     CSUBROUTINE---CHECK FOR ERRORS')                                           
      RETURN                                                                    
      END