
	program coord

c   ===================================================================

c  Yves Langevin, 1989
c		Updated StŽphane ERARD, juin 95

c	Compute coordinates (for second order even device) and angles.
c	Coordinates for other devices can be computed by changing the value
c	 of variable ordre.

c	Uses orbital parametres in file orbit.dat. Must be linked with 
c	 Orbit.for.

c	Longitudes are computed using the terrestrial convention 
c	 (from -180 to 180 eastward). This is required to write the data
c	 as 16-bits integers.
c	All angles are given in degrees (*100 in the output file).

c   ===================================================================

	implicit none

	integer*2 	tmois,i,id,imo,nspe,imois,ij,
     &			k,nbf,iph,ips,ipm,ipd,gain0,
     &			ih,im,is,jh,jm,js,jtic,
     &			Mirmin,x,y,xold,ymax,long,mir(300,26),
     &			ordre,buff(72),buff1(14)

	real*4		xp(6),yp(6),zp(6),ts,temps(300,26),
     &			xx,yy,zz,axi,exi,ainc,gomeg,pomeg,alo,ala,tobs,
     &			ys,xs,zs,xm,ym,zm,dpix,dpix1,dpix2,transm,
     &			rm,dani,ta,tb,am,an,an1,r,
     &			blong(4),blat(4),trans,dtmor,alongm,alatm,
     &			raym,durm(11),olong(4),olat(4),
     &			pinc,pemer,phase,ps,drm,dxm,dym,dzm
			

	character*55 	nom,ch2,ch
	character*15	fic
	character*25	nomf

	data durm/31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30./



558	trans=acos(-1.)/180.
	dtmor=18000.+34.*60+37
	call iniorb(0.)
50	type *,' month :'
	accept *,imois
	if(imois.eq.0)goto 558
	type *,' day :'
	accept *,ij
	ch='$1$DUA1:[PCCOMMON.BASEISM.CALIB]'
	open(unit=4,file=ch//'orbit.dat',status='old',shared)
	do 5120 i=1,12
	  read(4,5200)id,imo,ih,im,is,jh,jm,js,nom,nomf,nspe,gain0
	  read(4,5201)axi,exi,ainc,gomeg,pomeg,iph,ipm,ips,ipd,alo,ala
	  if(imois.eq.imo.and.ij.eq.id)goto 5121
5120	continue
5200	format(1x,2i3,2x,6i3,1x,a13,a13,i5,i2)
5201	format(1x,f8.2,f9.6,f9.4,2f8.3,4i3,2f4.1)
5121	close(unit=4)
	write(6,5200)id,imo,ih,im,is,jh,jm,js,nom,nomf,nspe,gain0
	write(6,5201)axi,exi,ainc,gomeg,pomeg,iph,ipm,ips,ipd,alo,ala
	ainc=ainc*trans
	gomeg=gomeg*trans
	pomeg=pomeg*trans	
	call ell(axi,exi,gomeg,pomeg,ainc)
	ts=0.
	if(imois.gt.1)then
	  do i=1,imois-1
	   ts=ts+durm(i)
	  enddo
	endif
	tmois=ts
	tobs=(tmois+ij-1)*86400.-dtmor
	ts=(((((ts+ipd-1.)*24.)+iph)*60.)+ipm)*60.+ips	
	type *,'   reference time for pericentre : ',ts
	call iniorb(ts)
	call posmar(ts,raym)
	type *,' Sun - mars distance:',raym
	call dir(0.,0.,xs,ys,zs)
	dpix=0.1*trans
	dpix1=1.1*dpix
	dpix2=0.92*dpix
	rm=3394.
	transm=trans*0.01308333

	ch2 ='$1$DUA1:[PCCOMMON.BASEISM.CAL_DATA.MARS]'

	nomf='/'
	type*,nom
	if(nom(1:9).eq.'mars0802.')then
		nomf='fond0802.dat'
		gain0=1
		ih=11
		im=05
		is=21
		jh=11
		jm=12
		js=08
		ymax=9
		xold=296
		mirmin=12
		raym=1.54679
		fic='pav'
	endif
	if(nom(1:9).eq.'mars1102.')then
		nomf='fond1102.dg2'
		gain0=2
		ih=16
		im=51
		is=43
		jh=16
		jm=58
		js=00
		ymax=9
		xold=274
		mirmin=233
		raym=1.55088
		fic='bib'
	endif
	if(nom(1:9).eq.'mars0103.')then
		nomf='fond0103.dg2'
		gain0=2
		ih=10
		im=45
		is=28
		jh=11
		jm=12
		js=36
		ymax=26
		xold=121
		mirmin=222
		raym=1.572447
		fic='syr'
	endif
	if(nom(1:9).eq.'mars1303.')then
		nomf='fond1303.dg2'
		gain0=2
		ih=11
		im=55
		is=25
		jh=12
		jm=15
		js=54
		ymax=26
		xold=91
		mirmin=156
		raym=1.587290
		fic='oly'
	endif
	if(nom(1:9).eq.'mars0703.')then
		nomf='fond0703.dat'
		gain0=1
		ih=3
		im=0
		is=47
		jh=3
		jm=27
		js=53
		ymax=26
		xold=121
		mirmin=2
		raym=1.579486
		fic='vmc'
	endif
	if(nom(1:9).eq.'mars2102.')then
		nomf='fond2102.dg2'
		gain0=2
		ih=10
		im=23
		is=3
		jh=10
		jm=44
		js=58
		ymax=26
		xold=98
		mirmin=228
		raym=1.562730
		fic='ara'
	endif
	if(nom(1:9).eq.'mars2702.')then
		nomf='fond2702.dg2'
		gain0=2
		ih=2
		im=41
		is=42
		jh=3
		jm=8
		js=51
		ymax=26
		xold=121
		mirmin=26
		raym=1.570081
		fic='dae'
	endif
	if(nom(1:9).eq.'mars1203.')then
		nomf='fond1403.dg2'
		gain0=2
		ih=3
		im=23
		is=42
		jh=3
		jm=50
		js=47
		ymax=26
		xold=121
		mirmin=2
		raym=1.585198
		fic='aur'
	endif
	if(nom(1:9).eq.'mars2603.')then
		nomf='fond2603.dg2'
		gain0=2
		ih=15
		im=03
		is=59
		jh=15
		jm=30
		js=52
		ymax=26
		xold=120
		mirmin=198
		raym=1.600988
		fic='heb'
	endif
	if(nom(1:9).eq.'mars2103.')then
		nomf='fond2103.dg3'
		gain0=3
		ih=12
		im=20
		is=36
		jh=12
		jm=47
		js=38
		ymax=26
		xold=120
		mirmin=162
		raym=1.595524
		fic='asc'
	endif
	if(nom(1:9).eq.'mars1403.')then
		nomf='fond1403.dg2'
		gain0=2
		ih=11
		im=50
		is=24
		jh=12
		jm=17
		js=25
		ymax=26
		xold=120
		mirmin=210
		raym=1.587860
		fic='gor'
	endif
	if(nom(1:9).eq.'phob2503.')then
		type*,'Mars files only '
		stop
	endif
	if(fic.eq.' ')then
		write(6,*)'Pb date'
		goto558
	endif
	type *,' Sun - mars distance:',raym

c	----- Reads time and mirror angle in the even channels file -----

	nbf=0
	long=36

	open(unit=8,file=ch2//fic//'even.cal',form='unformatted',
     *	recl=long,access='direct',status='old')

	do x=1,xold
	 do y=2,ymax

	  nbf=nbf+1
	  read(8'nbf)buff

	  jh=buff(1)
	  jm=buff(2)
	  js=buff(3)
	  jtic=buff(4)
	  temps(x,y)=jh*3600.+jm*60.+js+jtic*0.125
	  if (buff(5).ne.x)type*,'Lézard en x', x, y
	  if (buff(6).ne.y)type*,'Lézard en y', x, y
	  mir(x,y)=buff(8)

	 enddo
	enddo
	close(8)
	close(9)
	if(nbf.eq.0) type*, 'Timing problem...'

	type*,'End in line',xold,', sample',ymax
	type*,'Nb of pixels read:',nbf
	Type*,'Computing coordinates...'


c	**** Open output file ****

	long=7

	open(unit=2,file=fic//'_coor.dat',form='unformatted',
     *	recl=long,access='direct',status='new')

	nbf=0
	do x=1,xold
	 do y=2,ymax
	  nbf=nbf+1

	  if(temps(x,y).eq.0) then	!Blank record for missing pixels
	   pinc=0.
	   pemer=0.
	   phase=0.
	   mir(x,y)=0
	   do k=1,4
	    blat(k)=0.
	    blong(k)=0.
	   enddo
	   goto 174
	  endif


c	----- Computes the coordinates and angles -----


c		Correct viewing direction for the 4 devices

	 ordre=2	! Originaly computes for 2nd order even only

	  if (ordre.eq.4)dani=0.		!1st order odd
	  if (ordre.eq.2)dani=trans*0.20	!2nd order even
	  if (ordre.eq.1)dani=-trans*0.15	!2nd order odd
	  if (ordre.eq.3)dani=trans*0.34	!1st order even
	  ta=temps(x,y)+tobs
	  tb=ta-ts
	  am=(mir(x,y)-2048.)*transm-dpix1
	  an=dpix2
	  an1=an+dani
	  call dir(am,an1,xp(1),yp(1),zp(1))
	  an=-an
	  an1=an+dani
	  call dir(am,an1,xp(4),yp(4),zp(4))
	  am=am+dpix1+dpix1
	  call dir(am,an1,xp(3),yp(3),zp(3))
	  an=-an
	  an1=an+dani
	  call dir(am,an1,xp(2),yp(2),zp(2))
	  am=am-dpix1
	  call dir(am,dani,xp(5),yp(5),zp(5))
	  call pos(1,tb,xx,yy,zz,r)
	  do k=1,4
	    call inter(ta,xx,yy,zz,xp(k),yp(k),zp(k),
     &      xm,ym,zm,alatm,alongm)
	    if(xm*xm.lt.(1.e-4).and.ym*ym.lt.(1.e-4))then
	     pinc=0.
	     pemer=0.
	     phase=0.
	     mir(x,y)=0
	     do i=1,4
	      blat(i)=0.
	      blong(i)=0.
	     enddo
	     goto 174
	    endif
	    blat(k)=aLATm+ala	! Hand correction (global shift)
	    bLONG(K)=aLONGm+alo
	  enddo

	  call inter(ta,xx,yy,zz,xp(5),yp(5),zp(5),xm,ym,zm,alatm,alongm)
	  alongm=alongm+alo
	  alatm=alatm+ala
	  dxm=xx-xm		 		!Spacecraft-surface vector
	  dym=yy-ym
	  dzm=zz-zm
	  drm=sqrt(dxm*dxm+dym*dym+dzm*dzm)	!surface-spacecraft distance
	  pinc=(xs*xm+ys*ym+zs*zm)/rm
	  if(pinc.gt.(-1.))then
	    pinc=180.-acos(pinc)/trans		!Incidence
	  else
	    pinc=0.
	  endif
	  pemer=(dxm*xm+dym*ym+dzm*zm)/drm/rm
	  if(pemer.lt.1.)then
	    pemer=acos(pemer)/trans		!Emergence
	  else
	    pemer=0.
	  endif
	  phase=0.
	  ps=-(xs*dxm+ys*dym+zs*dzm)/drm
	  if(ps.lt.1.)then
	    phase=acos(ps)/trans		!Phase
	  else
	    phase=0.
	  endif

174	  if (mir(x,y-1).ne.0)then	!Stick the pixels together
	   blong(1)=olong(2)		! (cosmetic only)
	   blong(4)=olong(3)
	   blat(1)=olat(2)
	   blat(4)=olat(3)
	  endif
	  buff1(1)=x			! code in integer format
	  buff1(2)=y
	  do k=1,4
	   olong(k)=blong(k)
	   olat(k)=blat(k)
	   blong(k)=-blong(k)		!Terrestrial convention
	   if (blong(k).lt.-180.) blong(k)=360.+blong(k)
	   buff1(1+2*k)=int(blong(k)*100)
	   buff1(2+2*k)=int(blat(k)*100)
	  enddo
	  buff1(11)=int(pinc*100)
	  buff1(12)=int(pemer*100)
	  buff1(13)=int(phase*100)
	  buff1(14)=mir(x,y)
	  write(2'nbf)buff1


c	----- End of main loop -----

	 enddo
	enddo
	type*,'Nb of recorded pixels:',nbf
	type*,'The image is',xold,' lines x',ymax-1,' samples'
	close(2)

	end