SUBROUTINE ERRI3D(XI,YI,ZI,XJ,YJ,ZJ,VC1,DXLG,DYLG,DZLG,DXLA,DYLA,         
     1 DZLA,DIJ,AIJ,ZIJ,DAZ,XXII,ETA,P,L,VC2,J)                                 
C                                                                               
C        THIS SUBROUTINE WILL COMPUTE THE VARIANCE-COVARIANCE MATRIX OF THE     
C      SPATIAL DISTANCE AND ASTRONOMIC AZIMUTH AND ZENITH ANGLE AS COMPUTED     
C      IN THE INVERSE POSITIONING PROBLEM.                                      
C                                                                               
C        INPUT:                                                                 
C              XI,YI,ZI,XJ,YI,ZJ - CARTESIAN COORDINATES OF POINTS I AND J.     
C              VC1 - VARIANCE COVARIANCE MATRIX OF CARTESIAN COORDINATES        
C                    POINTS I AND J.                                            
C              DIJ,AIJ,ZIJ - COMPUTED INFORMATION - DIJ - SPATIAL DISTANCE.     
C                                                 - AIJ - ASTRONOMIC AZIMUTH.   
C                                                 - ZIJ - ZENITH ANGLE.         
C              DAZ - ASTRONOMIC AZIMUTH MINUS THE GEODETIC AZIMUTH.             
C              XXII,ETA - DEFLECTION OF THE VERTICAL COMPONENTS AT I.           
C              P,L - GEODETIC LATITUDE AND LONGITUDE OF POINT I.                
C        OUTPUT:                                                                
C              VC2 - VARIANCE COVARIANCE MATRIX OF THE SPATIAL DISTANCE,        
C                    ASTRONOMIC AZIMUTH AND ZENITH ANGLE.                       
C                                                                               
C                                           WRITTEN BY G.BOWIE,DEC.1977.        
C                                                                               
      IMPLICIT REAL*8(A-Z)                                                      
      DIMENSION J(3,6),JT(6,3),VC1(6,6),VC2(3,3),JVC(3,6)                       
      SP=DSIN(P)                                                                
      CP=DCOS(P)                                                                
      SL=DSIN(L)                                                                
      CL=DCOS(L)                                                                
      COTANA=DCOTAN(AIJ)                                                        
      J(1,1)=-(XJ-XI)/DIJ                                                       
      J(1,2)=-(YJ-YI)/DIJ                                                       
      J(1,3)=-(ZJ-ZI)/DIJ                                                       
      J(1,4)=-J(1,1)                                                            
      J(1,5)=-J(1,2)                                                            
      J(1,6)=-J(1,3)                                                            
      J(2,1)=DXLA*((DAZ*SP*CL+SL+ETA*CP*CL)- DTAN(AIJ)*(SP*CL-DAZ*SL+XXI        
     1I*CP*CL))/(DXLA**2+DYLA**2)                                               
      J(2,2)=DXLA*((DAZ*SP*SL-CL+ETA*CP*SL)-DTAN(AIJ)*(SP*SL+DAZ*CL+XXII        
     1 *CP*SL))/(DXLA**2+DYLA**2)                                               
      J(2,3)=DXLA*((-DAZ*CP-ETA*SL)-DTAN(AIJ)*(-CP+XXII*SP))/(DXLA**2+DY        
     1LA**2)                                                                    
      J(2,4)=-J(2,1)                                                            
      J(2,5)=-J(2,2)                                                            
      J(2,6)=-J(2,3)                                                            
      J(3,1)=-(XXII*SP*CL+ETA*SL-CP*CL+DCOS(ZIJ)*(XJ-XI)/DIJ)/DSQRT(DXLA        
     1 **2+DYLA**2)                                                             
      J(3,2)=-(XXII*SP*SL-ETA*CL-CP*SL+DCOS(ZIJ)*(YJ-YI)/DIJ)/DSQRT(DXLA        
     1 **2+DYLA**2)                                                             
      J(3,3)=-(-XXII*CP-SP+DCOS(ZIJ)*(ZJ-ZI)/DIJ)/DSQRT(DXLA**2+DYLA**2)        
      J(3,4)=-J(3,1)                                                            
      J(3,5)=-J(3,2)                                                            
      J(3,6)=-J(3,3)                                                            
      CALL TRNSD(JT,6,J,3,3,6)                                                  
      CALL MMULD(JVC,3,J,3,VC1,6,3,6,6)                                         
      CALL MMULD(VC2,3,JVC,3,JT,6,3,6,3)                                        
      RETURN                                                                    
      END