
	subroutine dir(am,an,dx,dy,dz)
c	 Direction d'observation dedimensionalisee (dx,dy,dz)
	common /ori/aix,aiy,aiz,ajx,ajy,ajz,akx,aky,akz
	cm=cos(am)
	sm=sin(am)
	cn=cos(an)
	sn=sin(an)
c	 (cm*cn,sn,cn*sm) est la direction de visee par rapport a la direction
c		antisolaire. On la ramene au repere stellaire
	dx=aix*cm*cn+ajx*sn+akx*cn*sm
	dy=aiy*cm*cn+ajy*sn+aky*cn*sm
	dz=aiz*cm*cn+ajz*sn+akz*cn*sm
	return
	end
	
	subroutine posmar(ts,raym)
c	 Distance Mars/Soleil (raym)
	pi=3.1415926535
	amar=1.5236915
	emar=0.09339494
	imar=1.8467848*trans
	gomar=49.160389*trans
	pomar=286.256111*trans
	tmar=141.25667*86400.
	permar=amar*sqrt(amar)*86400.*365.2425
	PARM=AMAR*(1.-EMAR*EMAR)
	c1=sqrt((1.-emar)/(1.+emar))
	c2=emar*sqrt(1.-emar*emar)
	dt0=(ts+tmar)/permar
	am0=dt0*2*pi
	am=am0
	dt=2.*atan(c1*tan(0.5*am))-c2*sin(am)/(1.+emar*cos(am))
	am=am+am0-dt	
	dt=2.*atan(c1*tan(0.5*am))-c2*sin(am)/(1.+emar*cos(am))
	am=am+am0-dt
	dt=2.*atan(c1*tan(0.5*am))-c2*sin(am)/(1.+emar*cos(am))
	am=am+am0-dt
	raym=parm/(1.+emar*cos(am))
	return
	end

	subroutine iniorb(ts)
c	 ts est le nb de s entre le 1/1/89 et le passage au periastre de 
c		l'orbite.
	common /ori/aix,aiy,aiz,ajx,ajy,ajz,akx,aky,akz
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	ror=9680.	
	por=28865.
	tm=88620.	!Periode rotation Mars
	rm=3394.	!rayon Mars
	rm2=rm*rm
	pi=3.1415926535
	trans=pi/180.
	perc=6001.5 !rm**1.5 * 2pi/racine(GM) pour calculer la per. de la sonde
c		Parametres orbitaux Mars
	amar=1.5236915		! en ua
	emar=0.09339494
	imar=1.8467848*trans	!en radians
	gomar=49.160389*trans
	pomar=286.256111*trans
	tmar=141.25667*86400.	!Nb de s entre pass au perihelie et le 1/1/89
	permar=amar*sqrt(amar)*86400.*365.2425	!Periode de revol en s
	PARM=AMAR*(1.-EMAR*EMAR)	!parametre Mars
	c1=sqrt((1.-emar)/(1.+emar))
	c2=emar*sqrt(1.-emar*emar)
	dt0=(ts+tmar)/permar
	am0=dt0*2*pi	!Anomalie moyenne de Mars
	am=am0
	dt=2.*atan(c1*tan(0.5*am))-c2*sin(am)/(1.+emar*cos(am))
	am=am+am0-dt	
	dt=2.*atan(c1*tan(0.5*am))-c2*sin(am)/(1.+emar*cos(am))
	am=am+am0-dt
	dt=2.*atan(c1*tan(0.5*am))-c2*sin(am)/(1.+emar*cos(am))
	am=am+am0-dt	!estimation iterative de l'anomalie vraie
	am0=am+pomar	!angle v/ noeud ascendant
	xox=0.95353810
	xoy=0.17452938
	xoz=0.24556993
	yox=-.08156049
	yoy=0.93421486
	yoz=-.34726141
	xcan=0.40464400		!Direction Canopus
	ycan=-.11682670	
	zcan=-.90698112
	cam0=cos(am0)
	sam0=sin(am0)
	aix=cam0*xox+sam0*yox
	aiy=cam0*xoy+sam0*yoy
	aiz=cam0*xoz+sam0*yoz
	delts=-aiz/trans
	alphs=atan(aiy/aix)
	if(aix.gt.0.)alphs=alphs+pi
	alphs=alphs/trans
	alx=-ycan*aiz+zcan*aiy
	aly=-zcan*aix+xcan*aiz
	alz=-xcan*aiy+ycan*aix
	am=sqrt(alx*alx+aly*aly+alz*alz)
	alx=alx/am
	aly=aly/am
	alz=alz/am
	amx=aiy*alz-aiz*aly
	amy=aiz*alx-aix*alz
	amz=aix*aly-aiy*alx
	cinc=cos(22.5*trans)
	sinc=sin(22.5*trans)
	akx=-amx*cinc-alx*sinc
	aky=-amy*cinc-aly*sinc
	akz=-amz*cinc-alz*sinc
	ajx=alx*cinc-amx*sinc
	ajy=aly*cinc-amy*sinc
	ajz=alz*cinc-amz*sinc
c	if(ts.gt.10.)then
c	  type *,' alphs,delts ',alphs,delts
c	  type *,' sonde '
c	  type *,aix,aiy,aiz
c	  type *,alx,aly,alz
c	  type *,amx,amy,amz
c	  type *,' instrument :'	!Orientation instrument v/ ref stellaire
c	  type *,aix,aiy,aiz	!au moment du periastre sonde
c	  type *,ajx,ajy,ajz
c	  type *,akx,aky,akz
c	endif
	return
	end

	subroutine pos(i,t,x,y,z,r)
c	 Vecteur sonde/centre de Mars (x,y,z,r)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	common /cell/ax,ex,par,perio,vamo,vaco,c1,c2,
     1  xa,ya,za,xb,yb,zb,bound(37)
	do k=1,37
	  if(bound(k).gt.t)goto 20
	end do
20	ang=10.*trans*(k-19+(t-bound(k-1))/(bound(k)-bound(k-1)))
	t1=time(ang)
910	r=par/(1.+ex*cos(ang))
	ang=ang+(t-t1)/r/r*vaco
	t1=time(ang)
	if(abs((t1-t)/t).gt.1.e-5)goto 910
	cang=cos(ang)
	sang=sin(ang)
	r=par/(1.+ex*cang)	!Dist sonde/ centre de Mars en rayons martiens
	x=r*rm*(cang*xa+sang*xb)	! Vecteur dimensionne
	y=r*rm*(cang*ya+sang*yb)
	z=r*rm*(cang*za+sang*zb)
	return
	end

	subroutine inter(t,x,y,z,dx,dy,dz,xm,ym,zm,alat,along)
c	 Coord du point observe (alat ,along)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	dist2=x*x+y*y+z*z-(x*dx+y*dy+z*dz)*(x*dx+y*dy+z*dz)
	if(dist2.ge.rm2)then
	  xm=0.
	  ym=0.
	  zm=0.
	  along=0.
	  alat=0.
	  return
	else
	  if(dist2.lt.0.)dist2=0.
	  ret=sqrt(rm*rm-dist2)
	  ps=-x*dx-y*dy-z*dz
	  xm=x+dx*(ps-ret)
	  ym=y+dy*(ps-ret)
	  zm=z+dz*(ps-ret)
	  alat=90.-acos(zm/rm)/trans
	  along=atan(ym/xm)/trans
	  if(xm.lt.0)along=along+180.   
	  along=wlong(t)-along
	  if(along.gt.360.)along=along-360.
	  if(along.lt.0.)along=along+360.
	  if(along.lt.0.)along=along+360.
	endif
	return
	end

	subroutine inter1(t,x,y,z,dx,dy,dz,xm,ym,zm,alat,along,dist,ang1,ang2)
c	 Coord du point observe (alat ,along), renvoie la dist au centre du
c		disque de la direction de visee.
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
c	 Distance mini de la direction de visee au centre Mars
	dist2=x*x+y*y+z*z-(x*dx+y*dy+z*dz)*(x*dx+y*dy+z*dz)
	dist=sqrt(dist2/rm2)	!dist direction d'obs-centre disque en rayons m.
	ang1=dist2/(x*x+y*y+z*z)	
	ang1=asin(sqrt(ang1))	!dist ang au centre
	ang2=atan((sqrt(dist2)-rm)/abs(x*dx+y*dy+z*dz))	!Dist ang au bord
	if(dist2.ge.rm2)then
	  xm=0.
	  ym=0.
	  zm=0.
	  along=0.
	  alat=0.
	  return
	else
	  if(dist2.lt.0.)dist2=0.
	  ret=sqrt(rm*rm-dist2)
	  ps=-x*dx-y*dy-z*dz
	  xm=x+dx*(ps-ret)
	  ym=y+dy*(ps-ret)
	  zm=z+dz*(ps-ret)
	  alat=90.-acos(zm/rm)/trans
	  along=atan(ym/xm)/trans
	  if(xm.lt.0)along=along+180.   
	  along=wlong(t)-along
	  if(along.gt.360.)along=along-360.
	  if(along.lt.0.)along=along+360.
	  if(along.lt.0.)along=along+360.
	endif
	return
	end

	function wlong(t)
c	Angle entre le meridien 0 et une direction cst (dir Soleil/aphelie M?)
c	t secondes apres le 01/01/89 a 0h (en degres)
	wlong=227.797835+350.891983*t/86400.
	iwl=(wlong+180.)/360.
	wlong=wlong-iwl*360.
	return
	end	  

	subroutine limbes(tp,am)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	common /cell/ax,ex,par,perio,vamo,vaco,c1,c2,
     1  xa,ya,za,xb,yb,zb,bound(37)
	call dir(am,0.,dx,dy,dz)
	dang=10.*trans
	ang=-pi-dang
	imod=1
	do k=1,37
	  ang=ang+dang
	  cang=cos(ang)
	  sang=sin(ang)
	  r=par/(1.+ex*cang)
	  x=r*(cang*xa+sang*xb)
	  y=r*(cang*ya+sang*yb)
	  z=r*(cang*za+sang*zb)
	  ps=x*dx+y*dy+z*dz
	  ds=x*x+y*y+z*z-ps*ps
	  if(k.eq.1)then
	    jmod=1
	    if(ds.gt.1.)jmod=2
	    imod=jmod
	  else
	    imod=1
	    if(ds.gt.1.)imod=2
	  endif	  
	  if(imod.ne.jmod)then
	    ang0=ang
	    es=ds
	    ang1=ang-(1.-ds)/(cs-ds)*dang
941	    cang=cos(ang1)	  
	    sang=sin(ang1)
	    r=par/(1.+ex*cang)
	    x=r*(cang*xa+sang*xb)
	    y=r*(cang*ya+sang*yb)
	    z=r*(cang*za+sang*zb)
	    ps=x*dx+y*dy+z*dz
	    fs=x*x+y*y+z*z-ps*ps
	    if(abs(fs-1.).gt.1.e-5)then
	      ang2=ang1
	      ang1=ang1-(1.-fs)/(es-fs)*(ang1-ang0)    
	      es=fs
	      ang0=ang2
	      goto 941
	    endif
	    ip=1
	    if(ps.lt.0.)ip=-1
	    t1=time(ang1)
	    type *,' time, ang, ip, imod : ',t1,(ang1/trans),ip,imod
	    xm=x-ps*dx
	    ym=y-ps*dy
	    zm=z-ps*dz
	    alat=90.-acos(zm)/trans
	    along=atan(ym/xm)/trans
	    if(xm.lt.0)along=along+180.   
	    along=wlong(tp+t1)-along
	    if(along.gt.360.)along=along-360.
	    if(along.lt.0.)along=along+360.
	    if(along.lt.0.)along=along+360.
	    type *,' pos sur mars : ',xm,ym,zm
	    type *,' long, lat    : ',along,alat
	  endif
	  cs=ds
	  jmod=imod
	end do  	 
	return
	end

	subroutine ell(axi,exi,gomeg,pomeg,ainc)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	common /cell/ax,ex,par,perio,vamo,vaco,c1,c2,
     1  xa,ya,za,xb,yb,zb,bound(37)
	xi=cos(gomeg)
	yi=sin(gomeg)
	zi=0.
	xj=-sin(gomeg)*cos(ainc)
	yj=cos(gomeg)*cos(ainc)
	zj=sin(ainc)
	cosp=cos(pomeg)
	sinp=sin(pomeg)
	xa=cosp*xi+sinp*xj
	ya=cosp*yi+sinp*yj
	za=sinp*zj
	xb=-sinp*xi+cosp*xj
	yb=-sinp*yi+cosp*yj
	zb=cosp*zj
	ra=rm+alt
	ax=axi
	ax1=ax/rm
	ex=exi
	c1=sqrt((1.-ex)/(1.+ex))
	c2=c1*(1.+ex)
	par=ax1*c2*c2
	perio=ax1*sqrt(ax1)*perc
	vaco=2.*pi*sqrt(par)/perc
	bound(1)=-0.5*perio
	bound(37)=0.5*perio
	dideg=10.*trans
	ang1=-18.*dideg
	do i=2,36
	  ang1=ang1+dideg
	  bound(i)=time(ang1)	!Intervalles de tp necessaires pour parcourir
	end do			!10 degres en anomalie vraie
	return
	end

	function time(ang)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	common /cell/ax,ex,par,perio,vamo,vaco,c1,c2,
     1  xa,ya,za,xb,yb,zb,bound(37)
	time=2.*atan(c1*tan(0.5*ang))-c2*ex*sin(ang)/(1.+ex*cos(ang))
	time=time*perio/2./pi
	return
	end


	function airmass(ta,am,x,y,z,xs,ys,zs)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	call dir(am,0.,xp,yp,zp)
	call inter(ta,x,y,z,xp,yp,zp,xm,ym,zm,alat,along)
	dxm=x-xm
	dym=y-ym
	dzm=z-zm
	drm=sqrt(dxm*dxm+dym*dym+dzm*dzm)
	psm=(xs*xm+ys*ym+zs*zm)/rm
	if(psm.gt.(-1.))then
	  ainc=-0.5*pi+acos(psm)	!angle v/ au plan horizontal
	else
	  ainc=0.5*pi
	endif
	psm=(dxm*xm+dym*ym+dzm*zm)/drm/rm
	if(psm.lt.1.)then
	  emer=0.5*pi-acos(psm)		!angle v/ au plan horizontal
	else
	  emer=0.5*pi
	endif
	sinai=sin(ainc)			!cos i
	sinem=sin(emer)			!cos e
	if(sinai.eq.-1..or.sinem.eq.0.)then
	 airmass=1.
	 return
	else
	 airmass=0.5*(1./sinai+1./sinem)
	endif
	return
	end

	function aincid(ta,am,x,y,z,xs,ys,zs)
	common /geo/pi,trans,por,ror,rm,rm2,tm,perc
	call dir(am,0.,xp,yp,zp)
	call inter(ta,x,y,z,xp,yp,zp,xm,ym,zm,alat,along)
	dxm=x-xm
	dym=y-ym
	dzm=z-zm
	drm=sqrt(dxm*dxm+dym*dym+dzm*dzm)
	psm=(xs*xm+ys*ym+zs*zm)/rm
	if(psm.gt.(-1.))then
	  ainc=-0.5*pi+acos(psm)
	else
	  ainc=0.5*pi
	endif
	psm=(dxm*xm+dym*ym+dzm*zm)/drm/rm
	if(psm.lt.1.)then
	  emer=0.5*pi-acos(psm)
	else
	  emer=0.5*pi
	endif
	aincid=sin(ainc)
	return
	end