c ----calcule coordonnees approchees + fonction horloge recepteur
c
c **********************************************************************
c ************************ initialisation ******************************
c **********************************************************************
c
c ----variables parametrees
c
	parameter(nd=200,npar=60)
	real*8 psd(nd),calcul(nd),epoch(nd)
	real*8 param(npar),mat(npar,npar),cppi(npar),x(npar)
	real g(nd,npar)
	integer*2 sat(nd),nupar(nd)
c
c ----parametres de navigation
c
	real*8 nav(32,28),orb(28),iono(8)
c
c ----chaines de caracteres et variables logiques
c
	character*2 typ(6)
	character*1 rep
	character*3 rcp
	character*12 fich
	character*80 ligne
	logical getsv(32),elimsv(32),elim
c
c ----autres variables
c
	integer kobs,kpar,stt(32),nb,nep
	integer ntyp,ityp,sv(10),nmes,sa,niter
	integer iyear,imonth,iday,ihour,imin
	real*8 xyz(3),week,sec,cel,t,pseud,t1,kr1(10),xlim
	real*8 xs(3),rel,calc,kr,ks,te,data(10),aa,res,cdd2,res1
c
c ----initialisation
c
c ----valeur approchee des coefficients d'ionosphere
c
c	iono(1)=5.d-9
c	iono(2)=0.d0
c	iono(3)=-4.d-8
c	iono(4)=0.d0
c	iono(5)=8.5d4
c	iono(6)=0.d0
c	iono(7)=-1.d5
c	iono(8)=0.d0
c
c ----ouverture du fichier de commande
c
	open(9,file='pseudo1.out')
c
	open(12,file='pseudo1.dat',status='old')
 510	read(12,'(a12)',end=511) fich
c
	cel=2.99792458d8
	do 101 i=1,32
	stt(i)=0
	elimsv(i)=.false.
 101	getsv(i)=.false.
c
c **********************************************************************
c ***************** ouverture des fichiers *****************************
c **********************************************************************
c
c ----fichier de sortie des resultats
c
	write(9,'(a,a)') ' Station : ',fich(1:4)
c
c ----fichier d'ephemerides et chargement des ephemerides
c
	open(6,file=fich(1:11)//'n',status='old')
	do 102 i=1,4
 102	read(6,'(a)') ligne
 103	continue
	read(6,'(a)',end=104) ligne
	read(ligne,'(i2)') sa
	if(.not.getsv(sa)) then
	getsv(sa)=.true.
	read(ligne,'(2x,5i3,f5.1,3d19.12)')
     +iyear,imonth,iday,ihour,imin,sec,nav(sa,2),nav(sa,3),nav(sa,4)
	call timconv(week,nav(sa,1),iyear,imonth,iday,ihour,imin,sec)
	do 105 i=1,6
 105	read(6,11) (nav(sa,4*i+l),l=1,4)
	else
	do 106 i=1,6
 106	read(6,'(a)') ligne
	endif
	goto 103
 104	close(6)
c
c ----fichier de donnees
c
	open(8,file=fich,status='old')
c
c **********************************************************************
c ********** premiere lecture du fichier pour information **************
c **********************************************************************
c
c ----lecture de l'entete d'un fichier RINEX
c
	do 201 i=1,6
 201	read(8,'(a)') ligne
c
c ----identifie le type de recepteur
c
	rcp=ligne(21:23)
	if(rcp.eq.'SER') then
	write(9,'(a)') 'Recepteur SERCEL NR52 ou TR5S'
	else if(rcp.eq.'ASH') then
	write(9,'(a)') 'Recepteur ASHTECH L-XII'
	else if(rcp.eq.'TRI') then
	write(9,'(a)') 'Recepteur TRIMBLE'
	else
	write(9,'(a)') 'Recepteur inconnu'
	endif
	read(8,'(a)') ligne
	read(8,'(3f14.4)') (xyz(i),i=1,3)
	read(8,'(a)') ligne
	read(8,'(a)') ligne
	read(8,'(a)') ligne
	read(ligne,'(i6)') ntyp
 	read(ligne,'(6x,6(4x,a2))') (typ(i),i=1,ntyp)
	do 202 i=1,ntyp
 202	if(typ(i).eq.'C1') ityp=i
c
c ----saut des lignes optionnelles
c
	do 203 i=1,22
 	read(8,'(a)') ligne
	if(ligne(61:61).eq.' ') goto 204
 203	continue
 204	nb=0
c
c ----satellites a eliminer pour le calcul
c
c	write(*,'(a)') ' Satellites a eliminer du calcul (0 = fin) : '
c 205	write(*,'(a,$)') ' Sv : '
c	read(*,*) j
c	if(j.ne.0) then
c	elimsv(j)=.true.
c	write(9,'(a,i2,a)') ' Satellite ',j,' elimine du calcul'
c	goto 205
c	endif
c
c ----lecture rapide des donnees pour comptage et statistique
c
	do 206 i=1,32
 206	getsv(i)=.false.
c
	do 207 nep=1,2000
	read(8,'(a)',end=208) ligne
	read(ligne,'(5i3,f11.7,3x,i3)')
     +iyear,imonth,iday,ihour,imin,sec,nmes
	read(ligne,'(32x,10i3)') (sv(i),i=1,nmes)
	do 207 i=1,nmes
	read(8,'(a)') ligne
	if(.not.elimsv(sv(i))) then
	getsv(sv(i))=.true.
	stt(sv(i))=stt(sv(i))+1
	endif
 207	continue
c
 208	write(9,'(a,i6)') ' Epoques de mesure : ',nep
	do 210 i=1,32
	if(getsv(i)) then
	write(9,'(a,2i5,a)') 'Satellite ',i,stt(i),' observations'
 	nb=nb+stt(i)
	endif
 210	continue
	write(9,'(a,i6)') ' Observations utilisables : ',nb
c
c **********************************************************************
c *********** determination des coordonnees approchees *****************
c **********************************************************************
c
c ----inialisation de variables pour l'inversion
c
	niter=0
	xlim=5.d5
	ipas=1+dble(nb)/dble(nd)
c
c ----rebobinage du fichier a chaque iteration
c
 215	continue
	elim=.false.
	k=0
	rewind(8)
 212	read(8,'(a)') ligne
	if(ligne(61:61).ne.' ') goto 212
c
c ----initialisation des boucles
c
	kobs=0
	kpar=3
	res=0.d0
	do 307 i=1,3
 307	param(i)=xyz(i)
c
c ----lecture d'un entete de bloc sur ipas
c
	do 301 j=1,nep
	read(8,'(a)',end=302) ligne
	read(ligne,'(5i3,f11.7,3x,i3)')
     +iyear,imonth,iday,ihour,imin,sec,nmes
	call timconv(week,t,iyear,imonth,iday,ihour,imin,sec)
	if(mod(j,ipas).eq.0) then
	read(ligne,'(32x,10i3)') (sv(i),i=1,nmes)
c
c ----lecture d'un bloc (boucle sur nmes mesures)
c
	kr=0.d0
	nelim=0
	do 303 i=1,nmes
	read(8,'(6(f14.3,2x))') (data(l),l=1,ntyp)
c
c ----test pour les donnees a eliminer
c
	if(elimsv(sv(i))) then
	elim=.true.
	nelim=nelim+1
	goto 303
	endif
c
c ----donnees a garder
c
	kobs=kobs+1
	psd(kobs)=data(ityp)
	epoch(kobs)=t
	sat(kobs)=sv(i)
c
c ----moyenne du decalage en temps de chaque observation
c
	do 304 l=1,28
 304	orb(l)=nav(sv(i),l)
	call direct1(orb,iono,xyz,t,xs,calc,rel)
	t=t-calc/cel
	call direct1(orb,iono,xyz,t,xs,calc,rel)
	te=dmod((t-orb(1)+3024.d2),6048.d2)-3024.d2
	ks=(orb(2)+orb(3)*te+orb(4)*te**2+rel-orb(27))*cel
	calcul(kobs)=calc-ks
	kr1(i)=psd(kobs)-calc+ks
	kr=kr+kr1(i)
c
c ----derivees partielles
c
	do 306 ii=1,3
 306	g(kobs,ii)=sngl((param(ii)-xs(ii))/calc)
c
c ----fin de lecture d'un bloc
c
 303	continue
c
c ----determination du kr moyen et rejet eventuel
c
	nmes=nmes-nelim
 	kr=kr/dble(nmes)
	res1=0.d0
	do 319 i=1,nmes
c	write(*,*) t,sv(i),kr1(i)
 319	res1=res1+(kr1(i)-kr)**2
	res1=dsqrt(res1/dble(nmes))
c
	if(res1.le.xlim) then
	kpar=kpar+1
	param(kpar)=kr
	do 305 l=kobs-nmes+1,kobs
 	calcul(l)=calcul(l)+kr
	res=res+(calcul(l)-psd(l))**2
 	g(l,kpar)=1.
 	nupar(l)=kpar
 305	continue
    	else
	write(*,*) 'elimine'
	kobs=kobs-nmes
	endif
c
c ----rejet de ipas-1 bloc sur ipas
c
	else
	do 321 i=1,nmes
 321 	read(8,'(a)') ligne
	endif
c
c ----fin de grande boucle de lecture des blocs
c
 301	continue
c
c **********************************************************************
c ********* debut de l'inversion pour coordonnees approchees ***********
c **********************************************************************
c
c ----ecarts type sur les donnees: 50m
c
 302	cdd2=50.d0**2
	res=dsqrt(res/dble(kobs))
	write(*,'(a,i6)') ' Nombre d'' observations : ',kobs
	write(*,'(a,i6)') ' Nombre de parametres   : ',kpar
	write(*,'(a,f9.2)') ' Residu moyen avant inversion : ',res
	res=0.d0
c
c ----remplissage de la seconde partie de MAT (termes lies aux
c		derivees partielles)
c
	do 308 i=1,kpar
	do 308 l=1,kpar
	aa=0.d0
	do 309 j=1,kobs
 	aa=aa+dble(g(j,i))*dble(g(j,l))
 309	continue
 	mat(i,l)=aa/cdd2
 308	continue
c
c ----remplissage de la premiere partie de MAT (incertitude sur les
c		parametres)
c
	do 317 i=1,kpar
 317	mat(i,i)=mat(i,i)+1.d0/(1.d4**2)
c
c ----calcul des variances a-posteriori sur les coordonnees
c
	do 310 j=1,3
	do 311 i=1,kpar
 311	x(i)=0.
	x(j)=1.
	call chol(mat,x,kpar)
 310	cppi(j)=dsqrt(dabs(x(j)))
c
c ----calcul du vecteur d'entree
c
 416	continue
c
	do 312 i=1,kpar
	x(i)=0.d0
	do 318 j=1,kobs
  	x(i)=x(i)+dble(g(j,i))*(psd(j)-calcul(j))/cdd2
 318	continue
 312	continue
c
c ----calcul du vecteur des variations
c
	call chol(mat,x,kpar)
c
c ----changement des parametres
c
	do 313 i=1,kpar
 313	param(i)=param(i)+x(i)
	do 326 i=1,3
 326	xyz(i)=param(i)
c
c ----residu moyen apres inversion
c
	res=0.d0
	do 314 j=1,kobs
	do 315 i=1,28
 315	orb(i)=nav(sat(j),i)
	t=epoch(j)
	call direct1(orb,iono,xyz,t,xs,calc,rel)
	t=t-calc/cel
	call direct1(orb,iono,xyz,t,xs,calc,rel)
	te=dmod((t-orb(1)+3024.d2),6048.d2)-3024.d2
	ks=(orb(2)+orb(3)*te+orb(4)*te**2+rel-orb(27))*cel
	kr=param(nupar(j))
	calcul(j)=calc-ks+kr
  	res=res+(psd(j)-calcul(j))**2
 314	continue
	res1=dsqrt(res/dble(kobs))
	write(*,'(a,f10.3,a)')
     +' Residu moyen apres inversion : ',res1,' m'
	res=dsqrt(x(1)**2+x(2)**2+x(3)**2)
	write(*,'(a,f10.3,a)') ' Ecart au point precedent : ',res,' m'
c
c ----recommence ?
c
c
	if(res.lt.2.d0) goto 415
c	write(*,'(a,$)') ' Autre iteration (o,n) ? '
c	read(*,'(a,$)') rep
c	if(rep.ne.'o') goto 415
 	niter=niter+1
	xlim=5.d0*res1
	if(niter.eq.1) goto 215
	if(.not.elim) then 
	kpar=3
	goto 416
	endif
	goto 215
c
c ----presentation des resultats a l'ecran
c
 415	write(*,*) 'Coordonnees et ecarts type:'
 	write(*,12)param(1),cppi(1)
	write(*,13)param(2),cppi(2)
	write(*,14)param(3),cppi(3)
c
c ----presentation des resultats dans le fichier trace
c
	write(9,'(/,a)') ' Coordonnees et ecarts type : '
	write(9,12) param(1),cppi(1)
	write(9,13) param(2),cppi(2)
	write(9,14) param(3),cppi(3)
	write(9,'(a,f10.2)') ' Residu moyen apres inversion : ',res1
	write(9,'(a,i6)') ' iterations : ',niter
c
c ----formats
c
 11	format(3x,4d19.12)
 12 	format(' x   : ',f13.2,' m     sigma x   : ',f11.2,' m')
 13 	format(' y   : ',f13.2,' m     sigma y   : ',f11.2,' m')
 14 	format(' z   : ',f13.2,' m     sigma z   : ',f11.2,' m')
c
c **********************************************************************
c ****** determination du decalage d'horloge a chaque epoque ***********
c **********************************************************************
c
c ----ouverture d'un fichier pour les decalages d'horloge recepteur
c
c	goto 504
	open(6,file=fich(1:11)//'h')
c
c ----rebobinage du fichier et saut d'eventuelles lignes optionnelles
c
	rewind(8)
 502	read(8,'(a)') ligne
	if(ligne(61:61).ne.' ') goto 502
c
c ----lecture des donnees du fichier Ashtech
c
 503	continue
c
c ----lecture d'entete de bloc
c
	read(8,'(a)',end=504) ligne
	read(ligne,'(5i3,f11.7,3x,i3)')
     +iyear,imonth,iday,ihour,imin,sec,nmes
	call timconv(week,t,iyear,imonth,iday,ihour,imin,sec)
	read(ligne,'(32x,10i3)') (sv(i),i=1,nmes)
	t1=t
c
c ----lecture de bloc
c ----calcul de la distance satellite - terre au temps de reception
c ----temps de vol
c ----calcul de la distance satellite - terre au temps d'emission
c ----decalage d'horloge recepteur
c
	kr=0.d0
	do 505 i=1,nmes
	read(8,'(6(f14.3,2x))') (data(j),j=1,ntyp)
	pseud=data(ityp)
	do 506 j=1,28
 506	orb(j)=nav(sv(i),j)
	call direct1(orb,iono,xyz,t,xs,calc,rel)
	t=t-calc/cel
	call direct1(orb,iono,xyz,t,xs,calc,rel)
	te=dmod((t-orb(1)+3024.d2),6048.d2)-3024.d2
	ks=(orb(2)+orb(3)*te+orb(4)*te**2+rel-orb(27))*cel
	kr1(i)=(pseud-calc+ks)
	kr=kr+kr1(i)
c
c ----ecriture des donnees
c
 505	continue
	kr=kr/dble(nmes)
	res=0.d0
	do 507 i=1,nmes
 507	res=res+(kr1(i)-kr)**2
	res=dsqrt(res/dble(nmes))
	write(6,'(3f10.0)') t1,kr,res
	goto 503
c
c ----retour en tete pour autre calcul de pseudo
c
c	close(6)
 504	close(8)
c	close(9)
	goto 510
c
 511	end
