SUBROUTINE VINDI(AE,F,XLAT1,XLONG1,FAZ,LINE,                              
     $                  LAT2,LON2,BA)                                           
      DOUBLE PRECISION A,ABS,AE,ATAN,B,BA,BE,C,CA,C2SM,CS,COS,CU1,CFAZ,         
     $                 DSIN,DCOS,DATAN,DSQRT,DABS,DSIG,DL,DLON,FAZ,             
     $                 LAT2,LINE,LON2,NEWSIG,TAN,TU1,TWOSM,XLAT1,               
     $                 XLONG1,U1,USQ,PI,F,FL,SIN,SQRT,SU1,SIG1,                 
     $                  SFAZ,SA,SIG,SS,ATAN2,DATAN2                             
      DATA PI/3.1415926535D0/                                                   
C  STATEMENT FUNCTIONS:                                                         
      SIN(A) = DSIN(A)                                                          
      COS(A) = DCOS(A)                                                          
      TAN(A) = DSIN(A) / DCOS(A)                                                
      SQRT(A) = DSQRT(A)                                                        
      ABS(A) = DABS(A)                                                          
      ATAN(A) = DATAN(A)                                                        
      ATAN2(A,B) = DATAN2(A,B)                                                  
C$                                                                              
C NAME        VINDI                                                             
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     DIRECT GEODETIC PROBLEM                                           
C                                                                               
C AUTHOR/DATE P.DELORME 780830 USING ALGORITHM OF T. VINCENTY                   
C                                                                               
C EXTERNALS   DSIN,DCOS,DSQRT,DABS,DATAN,DATN2                                  
C                                                                               
C CALLING     CALL VINDI(AE,F,XLAT1,XLONG1,FAZ,LINE,LAT2,LON2,BA)               
C                                                                               
C PARAMETERS  AE =  ELLIPSOID SEMI-MAJOR AXIS (M)                               
C             F = INVERSE OF ELLIPSOID FLATTENING                               
C             XLAT1 = INPUT LAT OF WESTERN POINT (DEG)                          
C             XLONG1 = INPUT LONG OF WESTERN POINT (DEG)                        
C             FAZ = FORWARD AZIMUTH(DEG)                                        
C             LINE = GEODESIC DISTANCE (M)                                      
C             LAT2 = OUTPUT LAT OF EASTERN POINT (DEG)                          
C             LON2 = OUTPUT LONG OF EASTERN POINT (DEG)                         
C             BA = BACKWARD AZIMUTH (DEG)                                       
C                                                                               
C LANGUAGE    FORTRAN                                                           
C$                                                                              
C                                                                               
C  MISCELLANEOUS VALUES:                                                        
C                                                                               
      FL = 1D0 / F                                                              
      BE = AE * (1D0 - FL)                                                      
C                                                                               
C  CALCULATE THE REDUCED LATITUDE AND ITS TRIG FUNCTIONS:                       
C                                                                               
      TU1 = (1D0 - FL) * TAN(XLAT1 * PI / 180D0)                                
      U1 = ATAN(TU1)                                                            
      SU1 = SIN(U1)                                                             
      CU1 = COS(U1)                                                             
C                                                                               
C  CALCULATE THE ANGULAR DISTANCE ON THE SPHERE FROM THE EQUATOR TO P1:         
C                                                                               
      CFAZ = COS(FAZ * PI / 180D0)                                              
      SIG1 = ATAN(TU1 / CFAZ)                                                   
C                                                                               
C  CALCULATE THE AZIMUTH OF THE GEODESIC AT THE EQUATOR:                        
C                                                                               
      SFAZ = SIN(FAZ * PI / 180D0)                                              
      SA = CU1 * SFAZ                                                           
      CA = SQRT(1D0 - (SA**2))                                                  
C                                                                               
C  CALCULATE U SQUARE:                                                          
C                                                                               
      USQ = (CA**2) * ((AE**2 - BE**2) / (BE**2))                               
C                                                                               
C  CALCULATE A:                                                                 
C                                                                               
      A = 1D0 + ((USQ / 256D0) * (64D0 + USQ * (-12D0 + 5D0 * USQ)))            
C                                                                               
C  CALCULATE B:                                                                 
C                                                                               
      B = (USQ / 512D0) * (128D0 + USQ * (-64D0 + 37D0 * USQ))                  
C                                                                               
C  SIGMA'S FIRST APPROXIMATION:                                                 
C                                                                               
      SIG = LINE / (BE * A )                                                    
C                                                                               
C                                                                               
C  ****************  BEGINNING OF ITERATION  *******************                
C                                                                               
C                                                                               
100   CONTINUE                                                                  
      TWOSM = (2D0 * SIG1) + SIG                                                
      C2SM = COS(TWOSM)                                                         
      CS = COS(SIG)                                                             
      SS = SIN(SIG)                                                             
C                                                                               
C  CALCULATE DELTA SIGMA:                                                       
C                                                                               
      DSIG = B * SS * (C2SM + 0.25D0 * B * CS * (-1D0 + 2D0                     
     $              * (C2SM**2)))                                               
      NEWSIG = (LINE / (BE * A)) + DSIG                                         
      IF(ABS(NEWSIG - SIG) .LE. 10D-10)  GOTO 70                                
      SIG = NEWSIG                                                              
      GOTO 100                                                                  
70    CONTINUE                                                                  
C                                                                               
C                                                                               
C  *******************  END OF ITERATION  ***************************           
C                                                                               
C                                                                               
C  CALCULATE THE LATITUDE OF EASTERN POINT:                                     
C                                                                               
      CS = COS(NEWSIG)                                                          
      SS = SIN(NEWSIG)                                                          
      LAT2 = ATAN2(((SU1 * CS) + (CU1 * SS * CFAZ)),((1D0 - FL)                 
     $            * SQRT((SA**2) + ((SU1 * SS) - (CU1 * CS * CFAZ))             
     $            **2))) * 180D0 / PI                                           
C                                                                               
C  CALCULATE THE DIFF IN LONGITUDE ON THE AUXILIARY SPHERE:                     
C                                                                               
      DL = ATAN2((SS * SFAZ),((CU1*CS) - (SU1 * SS * CFAZ)))                    
C                                                                               
C  CALCULATE C:                                                                 
C                                                                               
      C = (FL / 16D0) * (CA**2) * (4D0 + FL * (4D0 - 3D0 * (CA**2)))            
C                                                                               
C  CALCULATE THE DIFF IN LONGITUDE ON THE ELLIPSOID:                            
C                                                                               
      DLON = (DL - (1D0 - C) * FL * SA * (NEWSIG + C * SS * (C2SM               
     $        + C * CS * (-1D0 + 2D0 * (C2SM**2))))) * 180D0 / PI               
C                                                                               
C  CALCULATE THE LONGITUDE OF EASTERN POINT:                                    
C                                                                               
      LON2 = XLONG1 + DLON                                                      
C                                                                               
C                                                                               
C  CALCULATE THE BACKWARD AZIMUTH:                                              
C                                                                               
      BA = ATAN2((-1D0 * SA),((SU1 * SS) - (CU1 * CS * CFAZ)))                  
     $             * 180D0 / PI                                                 
      IF(BA .GE. 0D0)  GOTO 80                                                  
      BA = BA + 360D0                                                           
80    CONTINUE                                                                  
C                                                                               
C                                                                               
C                                                                               
C                                                                               
      RETURN                                                                    
      END