c ------------------------------------------------------
	subroutine direct1(orb,iono,xyzf,t,xs,calc,rel)
c ------------------------------------------------------
c
c ----calcul direct de distance terre-satellite + corrections
c
c ----variables utilisees par le sous-programme
c
	real*8 orb(28),iono(8),t,xs(3),calc,rel
	real*8 geogra(4),xyzf(3),ro,ionoc,elv,az,tropoc,cel
c
c ----initialisations
c
	cel=2.99792458d+8
c
c ----position du satellite a l'instant GPS t
c
	call possat(orb,xs,t,rel)
c
c ----distance … la station corrigee de la rotation de la terre
c
	call trajond(ro,xs,xyzf,cel)
c
c ----termes correctifs (ionosphere, troposphere)
c
	call angles(xs,xyzf,elv,az)
	call cargeo(xyzf,geogra)
	call coriono(iono,geogra,t,elv,az,ionoc)
	call cortropo(geogra,elv,tropoc)
c
c ----pseudo-range calcule sans les corrections d'horloge
c
	calc=ro+tropoc+ionoc
c
	end
c
c ------------------------------------
	subroutine possat(orb,x,t,rel)
c ------------------------------------
c
c ----position d'un satellite (donnee : t - retourne : x(3),rel)
c ----EN TEMPS GPS
c
	real*8 orb(28),x(3),t,rel,mu,ome,frel
	real*8 n0,tk,mk,ek,vk,phik,duk,dik,drk,uk
	real*8 rk,ik,xpk,ypk,omk,n,e,a
c
	frel=-4.442809305d-10
	mu=3.986005d14
	ome=7.2921151467d-5
c
	tk=t-orb(13)
	if(tk.gt.3024.d2) tk=tk-6048.d2
	if(tk.lt.-3024.d2) tk=tk+6048.d2
c
	e=orb(10)
	a=orb(12)**2
	n0=dsqrt(mu/(a**3))
	n=n0+orb(7)
	mk=orb(8)+n*tk
	rel=frel*dsqrt(a)*e*dsin(mk)
	ek=mk
	do 2 j=1,5
   2	ek=mk+e*dsin(ek)
	vk=2.d0*datan(dsqrt((1.d0+e)/(1.d0-e))*dtan(ek/2.d0))
	phik=vk+orb(19)
	ds2=dsin(2.d0*phik)
	dc2=dcos(2.d0*phik)
	duk=orb(9)*dc2+orb(11)*ds2
	drk=orb(18)*dc2+orb(6)*ds2
	dik=orb(14)*dc2+orb(16)*ds2
	uk=phik+duk
	rk=a*(1.d0-e*dcos(ek))+drk
	ik=orb(17)+dik+orb(21)*tk
	xpk=rk*dcos(uk)
	ypk=rk*dsin(uk)
	omk=orb(15)+(orb(20)-ome)*tk-ome*orb(13)
	x(1)=xpk*dcos(omk)-ypk*dcos(ik)*dsin(omk)
	x(2)=xpk*dsin(omk)+ypk*dcos(ik)*dcos(omk)
	x(3)=ypk*dsin(ik)
	end
c
c ---------------------------------
	subroutine trajond(r,x,y,c)
c ---------------------------------
c
c ----trajectoire de l'onde corrigee de la rotation de la terre
c
	real*8 x(3),y(3),w,c,r
	w=7.2921151467d-5
	r=dsqrt((x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2)
	r=r+(w/c)*(x(1)*y(2)-x(2)*y(1))
c
	end
c
c ---------------------------------------------
	subroutine cortropo(geogra,elev,tropoc)
c ---------------------------------------------
c
c ----modele standard de correction tropospherique (document OTAN)
c
	real*8 geogra(4),h,ns,dn,n1,c,dr,tropoc,fteta,elev
c
	h=geogra(2)/1.d3
	ns=3.248d2
	dn=-7.32d0*dexp(5.577d-3*ns)
	n1=ns+dn
	c=1.25d-1*dlog(n1/1.05d2)
	dr=7.32d-1+1.d-03*(n1/c)*(dexp(c*(1-h))-dexp(-8.d0*c))
	fteta=1/(dsin(elev)+0.00143/(dtan(elev)+0.0455))
	tropoc=dr*fteta
	end
c
c ------------------------------------------------------
	subroutine coriono(iono,geogra,temps,el,az,corr)
c ------------------------------------------------------
c
c ----modele standart de correction ionospherique (document OTAN)
c
	real*8 geogra(4),temps,el,az,corr,cel,pi
	real*8 lat,long,psi,phi,lambdai,phim,t,f,per,x,amp,tiono
	real*8 alpha(4),beta(4),iono(8)
c
c ----initialisations
c
	cel=2.99792458d+08
	pi=3.1415926535898d0
	do 3 i=1,4
   3	alpha(i)=iono(i)*1.d-9
	do 4 i=1,4
   4	beta(i)=iono(i+4)
c
	lat=geogra(4)/pi
	long=geogra(3)/pi
	psi=1.37d-2/((el/pi)+11.d-2)-22.d-3
	phi=lat+psi*dcos(az)
	if(lat.gt.416.d-3) phi=416.d-3
	if(lat.lt.-416.d-3) phi=-416.d-3
	lambdai=long+psi*dsin(az)/dcos(phi*pi)
	phim=phi+64.d-3*dcos((lambdai-1.617d0)*pi)
	t=4.32d+4*lambdai+temps+864.d2
	t=dmod(t,864.d2)
	f=1.d0+16.d0*((53.d-2-el/pi)**3)
	per=0.d0
	do 1 i=1,4
    1	per=per+beta(i)*(phim**(i-1))
	per=dmax1(per,72.d3)
	x=2.d0*pi*(t-504.d2)/per
	amp=0.d0
	do 2 i=1,4
    2	amp=amp+alpha(i)*(phim**(i-1))
	amp=dmax1(amp,0.d0)
	tiono=f*5.d-9
	if(dabs(x).lt.1.57d0) then
	tiono=tiono+f*amp*(1.d0-(x**2)/2.d0+(x**4)/24.d0)
	endif
	corr=tiono*cel
	end
c
c ----------------------------
	subroutine cargeo(x,r)
c ----------------------------
c
c ----passage des coordonnees cartesiennes aux geographiques
c ----x : tableau des coodonnees cartesiennes
c ----r : R,h,longitude,latitude
c
	real*8 x(3),r(4),a,e,p,n,pi
	pi=3.1415926535898d0
	a=6378135.d0
	e=0.0818191908d0
c
	r(1)=dsqrt(x(1)**2+x(2)**2+x(3)**2)
	r(3)=datan(x(2)/x(1))
	if(x(1).le.0) r(3)=r(3)+pi
	if(r(3).ge.pi) r(3)=r(3)-2.*pi
	p=dsqrt(x(1)**2+x(2)**2)
	r(4)=datan(x(3)/p/(1.d0-a*e*e/r(1)))
	do 1 i=1,5
	  n=a/dsqrt(1.d0-(e*dsin(r(4)))**2)
	  r(2)=p/dcos(r(4))-n
   1	  r(4)=datan(x(3)/p/(1.d0-n*e*e/(n+r(2))))
	end
c
c ----------------------------
	subroutine geocar(x,r)
c ----------------------------
c
c ----passage des coordonnees geographiques aux cartesiennes
c ----x : tableau des coordonnees cartesiennes
c ----r : R, h, longitude, latitude
c
	real*8 x(3),r(4),a,e,n
	a=6378135.d0
	e=0.0818191908d0
	n=a/dsqrt(1.d+00-(e*dsin(r(4)))**2)
	x(1)=(n+r(2))*dcos(r(4))*dcos(r(3))
	x(2)=(n+r(2))*dcos(r(4))*dsin(r(3))
	x(3)=(n*(1.d0-e*e)+r(2))*dsin(r(4))
	end
c
c ----------------------------------
	subroutine angles(s,g,el,az)
c ----------------------------------
c
c ----elevation et azimut du satellite(s) par rapport a la station(g)
c ----5/9/88

	real*8 s(3),g(3),el,az,pi,teta,rg,sg,psg,ps
	real*8 ua,ub,uc,va,vb,vc,wa,wb,wc,sga,sgb,sgc,psga,psgb,psgc
	pi=3.1415926535898d0
c
	rg=dsqrt(g(1)**2+g(2)**2+g(3)**2)
	wa=g(1)/rg
	wb=g(2)/rg
	wc=g(3)/rg
	va=-g(1)
	vb=-g(2)
	vc=-g(3)+rg/wc
	sg=dsqrt(va**2+vb**2+vc**2)
	va=va/sg
	vb=vb/sg
	vc=vc/sg
	ua=vb*wc-vc*wb
	ub=vc*wa-va*wc
	uc=va*wb-vb*wa
	sga=s(1)-g(1)
	sgb=s(2)-g(2)
	sgc=s(3)-g(3)
	sg=dsqrt(sga**2+sgb**2+sgc**2)
	ps=(wa*sga+wb*sgb+wc*sgc)/sg
	el=pi/2.-dacos(ps)
	psga=sga-sg*ps*wa
	psgb=sgb-sg*ps*wb
	psgc=sgc-sg*ps*wc
	psg=dsqrt(psga**2+psgb**2+psgc**2)
	psga=psga/psg
	psgb=psgb/psg
	psgc=psgc/psg
	teta=dacos(psga*ua+psgb*ub+psgc*uc)
	az=dacos(psga*va+psgb*vb+psgc*vc)
	if(teta.gt.(pi/2.)) az=2.*pi-az
c
	end
c
C ------------------------------------------------------------------
C     SOUS-PROGRAMME D'INVERSION DE MATRICE (CALCUL DE X AVEX Y=A*X)
C ------------------------------------------------------------------
C
      SUBROUTINE CHOL(MAT,X,NPAR)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT REAL*8 (M)
      PARAMETER (NFMX9=60)
      DIMENSION MAT(NFMX9,NFMX9),X(NFMX9),DIAG(NFMX9)
C
C ----------Inversion de la matrice MAT par la methode de CHOLEWSKI
C           Le vecteur de sortie est place dans X(I)
C
      DIAG(1)=DSQRT(MAT(1,1))
      DO 6 I=2,NPAR
    6 MAT(1,I)=MAT(I,1)/DIAG(1)
      DO 10 J=2,NPAR
      S=MAT(J,J)
      KMX=J-1
      DO 7 K=1,KMX
    7 S=S-MAT(K,J)**2
      DIAG(J)=DSQRT(S)
      IF(J.EQ.NPAR) GOTO 10
      IMI=J+1
      DO 8 I=IMI,NPAR
      S=MAT(I,J)
      DO 9 K=1,KMX
    9 S=S-MAT(K,I)*MAT(K,J)
      MAT(J,I)=S/DIAG(J)
    8 CONTINUE
   10 CONTINUE
      X(1)=X(1)/DIAG(1)
      DO 14 I=2,NPAR
      S=X(I)
      JMX=I-1
      DO 13 J=1,JMX
   13 S=S-MAT(J,I)*X(J)
   14 X(I)=S/DIAG(I)
      X(NPAR)=X(NPAR)/DIAG(NPAR)
      DO 16 K=2,NPAR
      I=NPAR+1-K
      JMI=I+1
      S=X(I)
      DO 15 J=JMI,NPAR
   15 S=S-MAT(I,J)*X(J)
   16 X(I)=S/DIAG(I)
c
      END
