SUBROUTINE ERRD3D(X,Y,Z,X1,Y1,Z1,P,L,R,DA,XI,ETA,AZ,ZE,VCDIR,VC,J)        
C                                                                               
C                                                                               
C        THIS ROUTINE COMPUTES THE VARIANCE-COVARIANCE MATRIX INVOLVED IN THE   
C      DIRECT PROBLEM OF GEODETIC POSITIONING IN A 3-DIMENSIONAL SYSTEM.        
C                                                                               
C        INPUT:                                                                 
C                                                                               
C              X,Y,Z - THE COORDINATES OF THE OBSERVATION STATION(IN METRES)    
C              X1,Y1,Z1 - THE COORDINATES OF THE OBSERVED STATION.(IN METRES)   
C              P,L - PHI,LAMDA OF THE OBSERVATION STATION IN THE GEODETIC       
C                    SYSTEM.(IN RADIANS)                                        
C              R - THE OBSERVED SPATIAL DISTANCE.(IN METRES)                    
C              AZ - THE OBSERVED AZIMUTH IN THE LOCAL ASTONOMIC SYSTEM.         
C                   (IN RADIANS)                                                
C              ZE - THE OBSERVED ZENITH ANGLE IN THE LOCAL ASTONOMIC SYSTEM.    
C                   (IN RADIANS)                                                
C              XI,ETA - THE DEFLECTION OF THE VERTICAL COMPONENTS AT THE        
C                       OBSERVATION STATION.(IN RADIANS)                        
C              VCDIR - THE VARIANCE-COVARIANCE MATRIX OF THE COORDINATES OF THE 
C                      OBSERVATION STATION AND OF THE ACTUAL OBSERVATIONS.      
C                                                                               
C        OUTPUT:                                                                
C                                                                               
C              VC - THE VARIANCE-COVARIANCE MATRIX OF THE COORDINATES OF THE    
C                   OBSERVATION STATION AND OF THE OBSERVED STATION.            
C                                                                               
C                                           WRITTEN BY G.BOWIE,DEC.1977.        
C                                                                               
      IMPLICIT REAL*8(A-Z)                                                      
      INTEGER I,K                                                               
      DIMENSION VC(6,6),VCDIR(6,6),JT(6,6),JVC(6,6),J(6,6)                      
      DO 1 I=1,6                                                                
      DO 1 K=1,6                                                                
      J(I,K)=0.D0                                                               
    1 CONTINUE                                                                  
      PI=3.141592653589793D0                                                    
      CP=DCOS(P)                                                                
      CL=DCOS(L)                                                                
      CA=DCOS(AZ)                                                               
      CZ=DCOS(ZE)                                                               
      SP=DSIN(P)                                                                
      SL=DSIN(L)                                                                
      SA=DSIN(AZ)                                                               
      SZ=DSIN(ZE)                                                               
      T1=-SZ*SA+DA*SZ*CA                                                        
      T2=DA*SZ*SA+SZ*CA                                                         
      T3=XI*SZ*SA-ETA*SZ*CA                                                     
      T4=CZ*CA+DA*CZ*SA-XI*SZ                                                   
      T5=-DA*CZ*CA+CZ*SA-ETA*SZ                                                 
      T6=-XI*CZ*CA-ETA*CZ*SA-SZ                                                 
      J(1,1)=1.D0                                                               
      J(2,2)=1.D0                                                               
      J(3,3)=1.D0                                                               
      J(4,1)=1.D0                                                               
      J(5,2)=1.D0                                                               
      J(6,3)=1.D0                                                               
      J(4,4)=(X1-X)/R                                                           
      J(4,5)=-R*(SP*CL*T1+SL*T2-CP*CL*T3)                                       
      J(4,6)=-R*(SP*CL*T4+SL*T5-CP*CL*T6)                                       
      J(5,4)=(Y1-Y)/R                                                           
      J(5,5)=-R*(SP*SL*T1-CL*T2-CP*SL*T3)                                       
      J(5,6)=-R*(SP*SL*T4-CL*T5-CP*SL*T6)                                       
      J(6,4)=(Z1-Z)/R                                                           
      J(6,5)=R*(CP*T1+SP*T3)                                                    
      J(6,6)=R*(CP*T4+SP*T6)                                                    
      CALL TRNSD(JT,6,J,6,6,6)                                                  
      CALL MMULD(JVC,6,J,6,VCDIR,6,6,6,6)                                       
      CALL MMULD(VC,6,JVC,6,JT,6,6,6,6)                                         
      RETURN                                                                    
      END