	program corpho

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

c	Yves Langevin, may 1990

c	Coregistrates spectral channels in the Phobos image-cube.
c	Results are interpolated raw data in DN.

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


	character*20	nomdec,nommtf
	dimension 	decal(128),pv(7)
	dimension 	x(24,25),y(24,25),z(128)
	integer*2 	jdat(24,25,128),kdat(24,25,128),
     &			m1(4),n1(4),ic(10),jc(10),icd(6),
     &			jcd(6),iv(7), nc0(6),nc1(6),nc2(6),
     &			ipixd(4),ipixf(4),spec(72),spec1(72),
     &			sto(24,25,8)

	character*55	ch, ch2

	data icd,jcd/-1,-1,0,1,1,0,0,-1,-1,0,1,1/
	data nc0,nc1,nc2/0,2,4,0,5,1,1,3,0,4,6,0,2,0,3,5,0,6/
	data ipixd/2,1,66,65/
        data ipixf/64,63,128,127/


c	Read scren coordinates
	ch2='$1$DUA1:[PCCOMMON.BASEISM.GEOMETRY]'
	open(unit=2,file=ch2//'pixloc.ctr',status='old')
	do i=1,24
	  do j=1,25
	    read(2,121)iscan,ielt,x(i,j),y(i,j)
	  end do
	end do
	close(unit=2)
121	format(1x,i2,i3,f8.2,f8.2)

c	Reads data corrected from telemetry errors
	ch='$1$DUA1:[PCCOMMON.BASEISM.EDT_DATA.PHOBOS]'
	open(1,file=ch//'phoodd.edt',form='unformatted',recl=36,
     &access='direct',status='old')
	open(2,file=ch//'phoeven.edt',form='unformatted',recl=36,
     &access='direct',status='old')
	nbf=0
	do i=1,24
	 do j=1,25
	  nbf=nbf+1
	  read(2'nbf)spec		!even channels
	  read(1'nbf)spec1		!odd channels

	  do k=1,8		!Saves instrumental parameters
	   sto(i,j,k)=spec(k)
	  enddo
	  do k=1,32
	   jdat(i,j,2*k)=spec(k+8)
	   jdat(i,j,2*k-1)=spec1(k+8)
	  enddo
	  do k=1,32
	   jdat(i,j,2*k+63)=spec(k+40)
	   jdat(i,j,2*k+64)=spec1(k+40)
	  enddo

	 enddo
	enddo
	type*,'nb of raw spectra read',nbf

	close(2)
	close(1)

716	format(3x,i3,3x,e11.5,3x,f6.4)
	do j=1,2
	  do i=1,24
	    do k=1,128
	      kdat(i,j,k)=jdat(i,1,k)
	    end do
	  end do
	end do

c	Read registration coefficients
	ch2='$1$DUA1:[PCCOMMON.BASEISM.CALIB]'
	nomdec='decal2.dat'
	open(unit=3,file=ch2//nomdec,status='old',shared)
	do k=1,128
	  read(3,716)i,decal(k),vv
	end do

	do i0=1,24
	do j0=3,25
	  x0=x(i0,j0)
	  y0=y(i0,j0)
	  i1=i0-3
	  if(i1.lt.1)i1=1
	  i2=i0+3
	  if(i2.gt.24)i2=24
	  j1=j0-3
	  if(j1.lt.1)j1=1
	  j2=j0+3
	  if(j2.gt.25)j2=25
	  do ip=1,128
	    x0k=x0+decal(ip)*0.6691	!shifts in the (u,v) plane of the image
	    y0k=y0+decal(ip)*0.7314
            dmin=3000.
 	    do i=i1,i2		!Looks for the 3 nearest neighbours
	      do j=j1,j2	! of the pixel's reference location 
	        dc=((x0k-x(i,j))*(x0k-x(i,j))
     &		   +(y0k-y(i,j))*(y0k-y(i,j)))
	        if (dc.lt.dmin)then
		  dmin=dc
		  ic(0)=i
		  jc(0)=j
	        endif
	      end do
	    end do
	    dx=x0k-x(ic(0),jc(0))
	    dy=y0k-y(ic(0),jc(0))

            do n=1,6		!Selects the 3 nearest neighbours defining
	      iv(n)=1		!a triangle that comprises the current pixel
	      ic(n)=ic(0)+icd(n)
	      jc(n)=jc(0)+jcd(n)
	      if(ic(n).lt.1.or.ic(n).gt.24)iv(n)=2
	      if(jc(n).lt.1.or.jc(n).gt.25)iv(n)=2
	      if(iv(n).eq.1)then
	        dxc=x(ic(n),jc(n))-x(ic(0),jc(0))
	        dyc=y(ic(n),jc(n))-y(ic(0),jc(0))
	        pv(n)=dxc*dy-dyc*dx
	      endif
	    end do
	    ic(7)=ic(1)
	    jc(7)=jc(1)
	    iv(7)=iv(1)
	    if(n.eq.1)pv(7)=pv(n)
	    amin=5000.
    	    do n=1,6
	      if(iv(n).eq.2.or.iv(n+1).eq.2)goto 1110
	      if((pv(n).ge.0.).and.(pv(n+1).le.0.))goto 1111
1110	    end do
	    if(iv(2).eq.1)then
	      n=1
	    else
	     n=4
	    endif
	    if(ic(0).eq.1.and.jc(0).eq.25)n=3
1111	    ic0=ic(nc0(n))
	    jc0=jc(nc0(n))
	    ic1=ic(nc1(n))
	    jc1=jc(nc1(n))
	    ic2=ic(nc2(n))
	    jc2=jc(nc2(n))
	  if(i0.eq.4.and.j0.eq.9)type *,ic0,jc0,ic1,jc1,ic2,jc2

c	Interpolate the channel's value at reference location
c	 from nearest neighbours
	    xc0=x(ic0,jc0)
	    yc0=y(ic0,jc0)
	    xc1=x(ic1,jc1)
	    yc1=y(ic1,jc1)
	    xc2=x(ic2,jc2)
	    yc2=y(ic2,jc2)
	    zc0=jdat(ic0,jc0,ip)+0.5
	    zc1=jdat(ic1,jc1,ip)+0.5
	    zc2=jdat(ic2,jc2,ip)+0.5

	    xy=xc0*yc2-xc2*yc0+xc2*yc1-xc1*yc2+xc1*yc0-xc0*yc1
	    yz=yc0*zc2-yc2*zc0+yc2*zc1-yc1*zc2+yc1*zc0-yc0*zc1
	    zx=zc0*xc2-zc2*xc0+zc2*xc1-zc1*xc2+zc1*xc0-zc0*xc1
	    xyz= (xc0*yc2-xc2*yc0)*zc1
     &	  	+(xc2*yc1-xc1*yc2)*zc0
     &		+(xc1*yc0-xc0*yc1)*zc2
	    kdat(i0,j0,ip)=(xyz-yz*x0k-zx*y0k)/xy
	  end do
	end do
	end do

c		----- Writes the corrected data -----

	nbf=0
	open(1,file='phoodd.res',form='unformatted',recl=36,
     &access='direct',status='new')
	open(2,file='phoeven.res',form='unformatted',recl=36,
     &access='direct',status='new')

	do i=1,24
	 do j=1,25
	  nbf=nbf+1

	  do k=1,8		!Saves instrumental parameters
	   spec(k)=sto(i,j,k)
	   spec1(k)=sto(i,j,k)
	  enddo
	  do k=1,32
	   spec(k+8)=kdat(i,j,2*k)
	   spec1(k+8)=kdat(i,j,2*k-1)
	  enddo
	  do k=1,32
	   spec(k+40)=kdat(i,j,2*k+63)
	   spec1(k+40)=kdat(i,j,2*k+64)
	  enddo

	  write(2'nbf)spec	!even channels
	  write(1'nbf)spec1	!odd channels

	 enddo
	enddo
	type*,'nb of registered spectra',nbf

	close(2)
	close(1)

 5000	stop
      	end

