	program photrace

c   ===================================================================
c               Yves Langevin, 1990
c               Updated july 1995, S. Erard

c       Plots spectral maps of Phobos on various devices.

c       Reads the data from the calibrated data files Phoeven/odd.cal
c        written by EtalPho.for (the name of the input directory is in ch2).

c       See file "SOFTINFO.TXT" for instructions. This program must be linked
c        with feneX.for, and possibly the Uniras library. It writes directly
c        in Regis language, so it can be used without Uniras on vt color
c        screens and regis/sixel printers (although lots of compilations
c        warnings would occur).


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


        implicit none


        integer*2       jpix1(10),jpix2(10),jpix(20),buff(72),buff1(72),
     &                  imois(10),ij(10),ivoi(5),buff2(14),long,hyp,
     &                  k,j,nbf,iph,ips,ipm,ipd,igain,ipx,k1,k2,
     &                  im,jh,jm,js,kh,km,ks,jtic,ipix1,ipix2,
     &                  xold,yold,id,nspe,jsup,ielim,iphob,gpcol,
     &                  ilis,ipuis,ipx1,ipx2,ititre,lsup,ksup,
     &                  ka,lk,ic1,ic2,ms,ns,itemp1,itemp2,
     &                  kv,ms1(5),ns1(5),j1,n1,iprod,xivmin

        integer*4       icol1(18),icol2(18),icol3(18),modcol(6),lun,
     &                  isort,ir,iolx,ioly,i,imode,ilun,ir1


        real*4          spec(128),blong(6),blat(6),x,y,dx,dy,
     &                  rapmin,rapmax,
     &                  matrice(30,26,128),angles(30,26,10),
     &                  temps(30,26),apx1,apx2,xmax,ymax,xmin,ymin,
     &                  ax,ex,ainc,gomeg,pomeg,alo,ala,tobs,
     &                  dpix,dpix1,dpix2,clong,clat,scirup,scirop,
     &                  xlim1,xlim2,xdeb,xran,xfin,ydeb,yran,yfin,
     &                  rm,dani,ta,tb,am,an,an1,r,airmass,ang,
     &                  dat1,dat2,scirap,eca,amoy,amoy2,blongm,blatm,
     &                  alongm,alatm,raym,pinc,pemer

        character*1     a,a1
        character*4     apix(10)

        character*55    nom,ch2,ch,ch1
        character*15    fic
        character*25    nomf
        character*14    nomfi
        character*72    tit,tiu,wdeb,wfin1,wfin2,wpol
        character*2     sh,sm,ss,th,tm,tt


	common /graf/xdeb,ydeb,xran,yran,xfin,yfin
	data icol1/ 360,360,340,320,300,270,240,220,205,
     &  180,172,165,157,150,142,135,127,120/
	data icol2/ 30,17*50 /
	data icol3/ 18*100/
	data modcol/2,2,2,2,2,2/


c       -----  Open and read input files  -----

        ch2='$1$DUA1:[PCCOMMON.BASEISM.CAL_DATA.PHOBOS]'
        ch='$1$DUA1:[PCCOMMON.BASEISM.GEOMETRY]'
        ch1='$1$DUA1:[PCCOMMON.BASEISM.CALIB]'

        long=36
        nbf=0
	fic='Pho'
	xold=24
	yold=26

        open(unit=8,file=ch2//fic//'even.cal',form='unformatted',
     *  recl=long,access='direct',status='old')
        open(unit=9,file=ch2//fic//'odd.cal',form='unformatted',
     *  recl=long,access='direct',status='old')
        open(unit=2,file=ch//fic//'_scrn.dat',form='unformatted',
     *  recl=7,access='direct',status='old')

        open(12,file='carte.dat',status='new')
        write(12,*)fic
        write(12,*)xold, yold

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

          nbf=nbf+1
          read(8'nbf)buff               !even channels
          read(9'nbf)buff1              !odd channels
          read(2'nbf)buff2              ! coordinates and angles

          jh=buff(1)
          jm=buff(2)
          js=buff(3)
          jtic=buff(4)
          temps(x,y)=jh*3600.+jm*60.+js+jtic*0.125      ! Time in s
          if (buff(5).ne.x)type*,'LŽzard en x', x, y
          if (buff(6).ne.y)type*,'LŽzard en y', x, y

          do i=1,32
           matrice(x,y,i*2)=buff(i+8)
           matrice(x,y,i*2-1)=buff1(i+8)
          enddo
          do i=1,32
           matrice(x,y,i*2+63)=buff(i+40)
           matrice(x,y,i*2+64)=buff1(i+40)
          enddo
          do i=1,128
           matrice(x,y,i)=matrice(x,y,i)/32767.0/2.
          enddo

          if (buff2(1).ne.x)type*,'LŽzard coord en x', x, y
          if (buff2(2).ne.y)type*,'LŽzard coord en y', x, y
          do k=1,10
           angles(x,y,k)=buff2(k+2)/100.
          enddo

         enddo
        enddo
        close(8)
        close(9)
        close(2)


c	raym=1.58
	call inigra

5002	ipx1=0
	type *,' Smoothing (1 = yes)?'
	accept *,ilis
936	type *,' First channel numbers:'
	accept 54,ipix1
	if(ipix1.le.0)goto 937
	ipx1=ipx1+1
c	dani=0.
c	if(ipix1.lt.65.and.mod(ipix1,2).eq.0)dani=trans*0.34
c	if(ipix1.gt.64.and.mod(ipix1,2).eq.1)dani=trans*0.34
	jpix1(ipx1)=ipix1
	goto 936
937	if(ipx1.eq.0)goto 5000
	apx1=1./float(ipx1)
	ipx2=0
938	type *,' Second channel numbers:'
	accept 54,ipix2
	if(ipix2.le.0)goto 939
	ipx2=ipx2+1
	jpix2(ipx2)=ipix2
	goto 938
54	format(i3)
945	format(i3,',')
946	format(i3,'/')
947	format(i3,' ')
939	ka=0
	if(ipx1.gt.1)then
	  do  k=1,(ipx1-1)
	    ka=ka+1
	    write(apix(ka),945)jpix1(k)
	  end do
	endif
	ka=ka+1
	write(apix(ka),946)jpix1(ipx1)
	if(ipx2.eq.0)goto 959
	apx2=1./float(ipx2)
	if(ipx2.gt.1)then
	  do k=1,(ipx2-1)
	    ka=ka+1
	    write(apix(ka),945)jpix2(k)
	  end do
	endif
	ka=ka+1
	write(apix(ka),947)jpix2(ipx2)
959	ka=ka+1
	do k=ka,10
	  apix(k)='    '
	end do
3500	xmax=150.
	xmin=240.
	ymin=-35.
	ymax=35.
	j1=jpix1(1)
	if(j1.lt.65)then
	  if(mod(j1,2).eq.0)then
	    dx=0.
	    dy=0.
	  else
	    dx=-4.5
	    dy=4.5
	  endif
	else
	  if(mod(j1,2).eq.0)then
	    dx=-2.5
	    dy=2.5
	  else
	    dx=2.
	    dy=-2.
	  endif
	endif
	imode=1
	type *,' Min value (<0 to change channels): '
	accept *, rapmin
	if(rapmin.lt.0.)goto 3000
	type *,' Maximum value: '
	accept *,rapmax
c	type *,' nivmin (180 default) :'
c	accept *,xivmin
c	xivmin=0.005
c	transm=trans*0.01308333
	itemp1=60*jh
	itemp1=60*itemp1+60*jm+js
	itemp2=60*kh
	itemp2=60*itemp2+60*km+ks
	type *,' titre ? (1 = yes)'
	accept *,ititre


c      --- Intialize graphics -----

	type *,'  vt340 screen (1), file (2), gpx (3), post (4) '
	accept *,isort
	call inifen(isort)
	xlim1=xdeb+1.5*xran
	xlim2=xdeb+2.*xran

	tit='    $'
	tiu(1:5)='    $'
	if(isort.gt.2)then
	  call gmslev('I','I','T','A')
	  CALL GOPEN
	  call gcharf('SOFT')
	  call gcharf('COMP')
	  if(modcol(isort).eq.1)then
	    call rcmode('BWS',16)
	    call ghalft
	    do 7445 i=2,16
7445	    call gcolor(i,(18-i),0,0,1)
	  else
	    call rcmode('HLS',100)
	    call gcolor(2,icol1,icol2,icol3,18)
	  endif
	else
	  open(unit=2,file=ch1//'tabcol.dat',status='old')
	  LUN=6
	  if(isort.eq.2)then
	    type *,' Background: white (1) black (2) ?'
	    accept *,ilun
	    LUN=LUN+ilun
	    type *,'  title : '
	    accept 2226,tit
2226	    format(a72)
	    open(unit=LUN,file='regis.dat',status='new')
	  endif
	  ir1=mod(LUN,2)
	  do lk=1,29
	    read(2,881)wdeb
	    write(LUN,881)wdeb
	  end do
	  if(LUN.eq.7)then
	    write(LUN,923)
923	    format(' S(I7,E)')
	  endif
	  read(2,881)wfin1
	  read(2,881)wfin2
	  iolx=0
	  ioly=0
881	format(a72)
	endif
	if(isort.gt.2)then
	 call axes(tit,tiu,xmin,xmax,ymin,ymax)
	else
	 call mirl(LUN)
	 call axesl(isort,LUN,xmin,xmax,ymin,ymax)
	endif	
8000	ic1=0
	ic2=0
	amoy=0.
	amoy2=0.
8500	format(3i6,f11.2)


c      ----  Map data -----

	do ms=1,24
	 do ns=1,25
c          if (angles(ms,ns,11).lt.0.) goto 6000
          ic1=ic1+1
          pemer=angles(ms,ns,10)*3.14159/180.   !Take the cosine
          pinc=angles(ms,ns,9)*3.14159/180.
          pinc=acos(pinc)
          pemer=acos(pemer)

	  dat1=1.
	  do k=1,ipx1
	    k1=jpix1(k)
	    dat1=dat1*matrice(ms,ns,k1)
	  end do

	  dat2=1.
	  if(ipx2.ne.0)then
	    do k=1,ipx2
	      k2=jpix2(k)
	      dat2=dat2*matrice(ms,ns,k2)
	    end do
	  endif

	  if(dat2.le.0..or.dat1.le.0.)then
	    scirap=0.001
	  else	
	    scirap=dat1**apx1/(dat2**apx2)
	  endif
	  if(scirap.lt.0.005)goto 6000
	  if(matrice(ms,ns,12).lt.0.005)goto6000
	  if(ilis.eq.1)then
	    ms1(1)=ms
	    ns1(1)=ns-1
	    ms1(2)=ms
	    ns1(2)=ns+1
	    ms1(3)=ms-1
	    ns1(3)=ns
	    ms1(4)=ms+1
	    ns1(4)=ns
	    scirup=2.*scirap
	    do kv=1,4
	      iprod=(ms1(kv)*ns1(kv)*(25-ms1(kv))*(26-ns1(kv)))
	      if(iprod.eq.0)then
	        scirup=scirup+scirap
	      else
	        dat1=1.
	        do k=1,ipx1
	          k1=jpix1(k)
	          dat1=dat1*matrice(ms1(kv),ns1(kv),k1)
	        end do
	        dat2=1.
	        if(ipx2.ne.0)then
	          do k=1,ipx2
	            k2=jpix2(k)
	            dat2=dat2*matrice(ms1(kv),ns1(kv),k2)
	          end do
	        endif
	        if(dat2.le.0..or.dat1.le.0.003)then
	          scirop=scirap
	        else	
	          scirop=dat1**apx1/(dat2**apx2)
	        endif
	        scirup=scirup+scirop
	      endif
	    end do
	    scirap=scirup/6.
	  endif

	  do 800 k=1,4
            if(angles(ms,ns,(k-1)*2+1).lt.10.)goto 6000
	    k1=k
	    n1=ns
	    if(angles(ms,ns-1,(k-1)*2+1).gt.10.)then
	      if(k.eq.1)then
	        n1=ns-1
	        k1=2
	      endif
	      if(k.eq.4)then
	        n1=ns-1
	        k1=3
	      endif
	    endif
	    clong=angles(ms,n1,(k1-1)*2+1)+dx
	    clat=angles(ms,n1,k1*2)+dy
	    call recadr(imode,clong,clat,
     &      xmin,xmax,ymin,ymax,blong(k),blat(k))	
800	  continue

	  if(modcol(isort).eq.1)then
	    ir=(scirap-rapmin)/(rapmax-rapmin)*14.9999+2
	    if(ir.lt.2)ir=2
	    if(ir.gt.16)ir=16
	  else
	    ir=(scirap-rapmin)/(rapmax-rapmin)*17.9999+2
	    if(ir.lt.2)ir=2
	    if(ir.gt.19)ir=19
	  endif
	  if(blong(3).gt.xlim1.and.blong(3).lt.xlim2)goto 666
	  if(blong(1).gt.xlim1)goto 6000
	  if(blong(3).gt.xfin)goto 6000
	  if(blong(1).lt.xdeb)goto 6000
	  if(blat(2).lt.ydeb)goto 6000
	  if(blat(4).gt.yfin)goto 6000
	  amoy=amoy+scirap
	  amoy2=amoy2+scirap*scirap
	  ic2=ic2+1
	  if(ir.ne.ir1)then
	    if(isort.gt.2)then
	      call rsurf(blong,blat,4,ir,0.0)
	    else
	      call lsurf(LUN,blong,blat,ir,iolx,ioly)
	    endif
	  endif
6000	 enddo
6010	enddo
666	if(ic2.ne.0)then
	  amoy=amoy/ic2
	  amoy2=amoy2/ic2
	endif
	eca=sqrt(abs(amoy2-amoy**2))
        write(sh,fmt=2221)jh
        write(sm,fmt=2221)jm
        write(ss,fmt=2221)js
        write(th,fmt=2221)kh
        write(tm,fmt=2221)km
        write(tt,fmt=2221)ks
2221	format(i2.2)
521	continue
	write(unit=tiu,fmt=2223)amoy,eca,(apix(k),k=1,9)
2223	format(f8.4,' +/-',f8.5,'  chan: ' 9a4,'$')
	call echell(isort,LUN,rapmin,rapmax)
	if(ititre.eq.1)call titre(isort,LUN,tiu,tit)
	call gricor(isort,LUN,xmin,xmax,ymin,ymax)
	if(isort.lt.3)then
	  if(LUN.ne.7)then 
            write(LUN,910)
910	  format(' W(I7)S(I0)')
	  else
	    write(LUN,911)
911	  format(' W(I0)S(I7)')
	  endif
	  write(LUN,912)
	  write(LUN,913)
	  write(LUN,914)
	  write(LUN,915)
	endif
912	format(' P[410,168]T(A0) ','''','+30','''')	
913	format(' P[460,415]T(A0) ','''','180','''')	
914	format(' P[343,415]T(A0) ','''','210','''')	
915	format(' P[194,415]T(A0) ','''','240','''')	
	if(isort.gt.2)then
	  CALL GCLOSE
	  accept 2441,a1
2441	  format(a1)
	else
	  close(unit=2)
	  write(LUN,881)wfin1
	accept 2441,a1
	  write(LUN,881)wfin2
	  if(LUN.eq.4)close(unit=4)
	endif
	type 2225,amoy,eca
2225	format(3x,' Mean: ',f9.5,' deviation: ',f9.5)
        write(12,*)tiu
        write(12,*)tit
	goto 3500

3000	close(12)
	type *,' Other channels? (y = yes)'
	accept 141,a
	if(a.eq.'y') then
	  go to 5002
	endif
141	format(a1)

5000	stop
      	end
