c program to fit photometric constants, modified May 22, 1997
c
c  major overhaul, Jun 18, 1997 to include plotting
c
c  modified Nov 12 1997 to allow using a reference filter for color
c    term which has no observed data (but have standard data)
c   the 8 filters now have fixed color based on the standard table
c     and are referred to by its number in the standard table
c
c   variables:         amt(8): default airmass terms 
c                      ct(8) : default color terms
c            common /ststar/ standard star table
c                      stnam: ascii name
c                      rmag(i,j): jth filter mag for ith star
c                      color(j):  ascii name for filter name
c                      alph(i):   right ascension in h.hhh 
c                      dec(i):    declination     in d.ddd
c                      nos=number of stars, noc=number of filters
c            common /stmag/: observed standard data
c                      stn(i): ascii name
c                      stmag(i,j): instr mag for jth filter for ith star
c                      idst(i): standard star table id number
c                      idco(j): color id
c                      ha(i)  : hour angle
c                      dec(i) : dec
c                      ts(j,i): integration time
c                      airmass(i): air mass
c                      nsd:             ncd
c                      iday, mon, iyr: day month year-1900
c                      obse: observatory name
        DIMENSION AMT(8),CT(8)
        CHARACTER*12 STNAM,STN
        CHARACTER*60 COMM,splot*18,hplot*18
        CHARACTER COLOR*1,obse*4,mnd*4,mtemp*4
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *                  ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *       TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
	common/plotdata/xx(200,8),yy(200,8),nps(8),a0(8),a1(8),a2(8),
     *                  itype,numo(200,8)
        common/plotdev/hplot,splot,ipinit
        write(6,*) ' Photometric Reduction Program: '
        DATA MSTOP/'q'/,MRSTF/'rstf'/,MSTDI/'istd'/,MARED/'ared'/
     *     , MTSTF/'tstf'/,MTDAT/'tdat'/,MSTDF/'rstd'/,MSTIME/'stim'/
     *     ,M2RED/'2red'/,MCRED/'cred'/,MPCRED/'pcre'/,MPARED/'pare'/
     *     ,MP2RED/'p2re'/,MPDAT/'pdat'/,MPSTF/'pstf'/,M0RED/'0red'/
     *     ,MP0RED/'p0re'/
        CALL RSTFIL(amt,ct,*99)
  99    write(6,9)
9       FORMAT( ' PRP>> ',$)
        NFP=0
        NFP2=0
        read(5,300) COMM
 300    FORMAT(A60)
        DO 310 ICH=1,10
          READ(COMM(ICH:ICH),301,END=311) MTEMP
          IF (MTEMP.EQ.' ') GO TO 311
 310    CONTINUE
 311    IF (COMM(1:1).EQ.'%') THEN
c          CALL LIB$SPAWN(COMM(2:60),,,,,'PHORED')
          GO TO 99
        END IF
	if (ich.le.1) goto 99
        READ(COMM(1:ICH-1),301) MND
 301    FORMAT(A4)
        READ(COMM(ICH:10),302,END=202,ERR=202) NFP,BLANK1,NFP2,BLANK2
 302    FORMAT(1X,I1,A1,I1,A1)
 202    IF (MND.EQ.'quit') GO TO 1000
        IF (MND.EQ.'help') GO TO 400
        IF (MND.EQ.'rstf') CALL RSTFIL(amt,ct,*99)
        IF (MND.EQ.'istd') CALL STANIN(*99,1,amt)
        IF (MND.EQ.'rstd') CALL STANIN(*99,0,amt)
        IF (MND.EQ.'0red') CALL REDUCE(*99,0,0,AMT,CT )
        IF (MND.EQ.'cred') CALL REDUCE(*99,1,0,AMT,CT )
        IF (MND.EQ.'ared') CALL REDUCE(*99,2,0,AMT,CT )
        IF (MND.EQ.'2red') CALL REDUCE(*99,3,0,AMT,CT )
        IF (MND.EQ.'p0re') CALL REDUCE(*99,0,1,AMT,CT )
        IF (MND.EQ.'pcre') CALL REDUCE(*99,1,1,AMT,CT )
        IF (MND.EQ.'pare') CALL REDUCE(*99,2,1,AMT,CT )
        IF (MND.EQ.'p2re') CALL REDUCE(*99,3,1,AMT,CT )
        IF (MND.EQ.'tstf') CALL TSTF(*99,0)
        IF (MND.EQ.'pstf') CALL TSTF(*99,1)
        IF (MND.EQ.'tdat') CALL TDAT(*99,0)
        IF (MND.EQ.'pdat') CALL TDAT(*99,1)
	if (mnd.eq.'plot') call plotdel(*99)
	if (mnd.eq.'calc') call calcat(*99)
	if (mnd.eq.'ph2s') call pho2stf(*99)
        IF (MND.EQ.'stim') CALL STIME(*99)
 400    write(6,*) ' AVAILABLE COMMANDS: quit,rstf,rstd,istd.tstf,'
     *          ,'pstf,tdat,pdat,stim,'
        write(6,*) '  2redu, aredu, credu, p2redu, paredu, pcredu'
     *                ,' 0redu,p0redu,plot,calcat,ph2s'
        GO TO 99
 1000   STOP
        END
C///////////////////////////////////////////////////////////////////////
C
c subroutine to read standard star file
C
        SUBROUTINE RSTFIL(amt,ct,*)
	dimension amt(8),ct(8)
        CHARACTER STNAM*12,FILEN*18,COLOR*1
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *           ,DEC(700),NOS,NCO
 601    write(6,1) 
 1      format(' Enter standard star catalog file name (n for none): '
     *         ,$)
        read(5,121,end=999) FILEN
	if (filen.eq.'n') return1
 121    FORMAT(A18)
        OPEN(UNIT=28,FILE=FILEN,status='OLD',ERR=602)
	goto 104
 602	write(6,*) '  *** no such file ***'
	goto 601
c read color
 104    READ(28,*,END=105) NCO,(COLOR(I),I=1,NCO)
 105    NOS=0
c read standard table data
        DO 101 I=1,700
          READ(28,*,END=100) STNAM(I), (RMAG(I,K),K=1,NCO),
     #                 IRH,IRM,RS,IDED,IDM,DS
          ALPH(I)=IRH*3600.+IRM*60.+RS
          IF (IDED.LT.0) then
	    IDM=-IDM
	    DS=-DS
	  else if (idm.lt.0) then
	    DS=-DS
	  end if
          DEC(I)=IDED*3600.+IDM*60.+DS
          NOS=NOS+1
 101	continue
c set up default colors and airmass
	do 102 i=1,nco
	  if (color(i).eq.'I'.or.color(i).eq.'i') then
	    amt(i)=-.037
	    ct(i)=-.03
	  else if (color(i).eq.'R'.or.color(i).eq.'r') then
	    amt(i)=-.052
	    ct(i)=.15
	  else if (color(i).eq.'V') then
	    amt(i)=-.11 
	    ct(i)=0
	  else if (color(i).eq.'g') then
	    amt(i)=-.11 
	    ct(i)=-.375
	  else if (color(i).eq.'B'.or.color(i).eq.'b') then
	    amt(i)=-.18  
	    ct(i)=-.02
	  else if (color(i).eq.'U'.or.color(i).eq.'u') then
	    amt(i)=-.300
	    ct(i)=-.09
	  else 
	    amt(i)=0.
	    ct(i)=.0
	  end if
 102	continue
 100    CLOSE(UNIT=28)
 999 	RETURN1
        END
c//////////////////////////////////
c  subroutine to input observatory data
c
	subroutine observatory(obse,rlong,rlat,utmzt)
	character obse*4
C  SET LONGITUDE(W), LATITUDE, AND UT-ST FOR OBSERVATORY
	if (obse.eq.' ') then
 600      write(6,1) 
 1	format(' Enter observatory name [kpno,ctio,cfht,lco,int]: ',$)
          read(5,'(a4)',END=999) OBSE
	end if
        IF (OBSE.EQ.'kpno'.OR.OBSE.EQ.'sokp') THEN
          RLAT=31.958
          RLONG=111.595
          UTMZT=7.
        ELSE IF (OBSE.EQ.'ctio') THEN
            RLAT=-30.1689
            RLONG=70.815
            UTMZT=4
        ELSE IF (OBSE.EQ.'cfht') THEN
              RLAT=19.8261
              RLONG=155.4717
        ELSE IF (OBSE.EQ.'int') THEN
            RLAT=28.75
            RLONG=17.8667
        ELSE IF (OBSE.EQ.'lco') THEN
            RLAT=-29.0083
            RLONG=70.7000
        ELSE
 6038     write(6,2)  rlon,rlat
 2	format('  No recognisable name supplied...'/
     *     '    Enter longitude and latitude of observatory [',2f7.1
     *        ,']: ',$)
	  call read2r(rlon,rlat,*999,*6038)
        END IF
 999	return
	end
C//////////////////////////////////////////////////////////////
c
C     SUBROUTINE TO READ STANDARD DATA
C             IPAR=0     READ FROM A FILE
C             IPAR=1     READ FROM TERMINAL INPUT
C
        SUBROUTINE STANIN(*,IPAR,amt)
        DIMENSION IOH(200),IOM(200),OS(200),amt(8)
        CHARACTER STNAM*12,STN*12,FILEN*18,ncdstr*1,fstr140*80,
     *             fstr141*80,fstr143*80,tstr*60
        CHARACTER COLOR*1,OBSE*4,comstr*60
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *           ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *        TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
c read header
        IF (IPAR.EQ.0) THEN   
	  obse=' '
c read from existing files
 660      write(6,3) 
 3	format('  Enter standard star data file name: ',$)
          read(5,11,END=999) FILEN
 11     FORMAT(A18)
          OPEN(UNIT=28,FILE=FILEN,status='OLD',ERR=660)
 502	READ(28,*,err=500) NCD,(IDCO(I),I=1,NCD),IDAY,MON,IYR
    	  call setairmass(amt,idco,color,ncd)
	  go to 501
c  read comment if first line is not date
 500	  backspace 28
	  read(28,19) comstr
 19	format(a60)
	  write(6,19) comstr
	  go to 502
 501	continue
          IF (IDAY.EQ.0) THEN
            IUTAIR=1
          ELSE
            IUTAIR=0
          END IF
	  call observatory(obse,rlong,rlat,utmzt)
          GO TO 200
        ELSE    
c input data from terminal
	  call observatory(obse,rlong,rlat,utmzt)
          write(6,4) 
 4	format(
     * ' Do you want to enter AIRMASS (a) or TIME OF OBS (t)? :',$)
          read(5,11,END=999) FILEN
          IF (FILEN.EQ.'t') THEN
 600        write(6,*) ' ENTER DATE (DD MM (YR-1900)) (UT) OF OBS'
            read(5,*,ERR=600,END=999) IDAY,MON,IYR
            IUTAIR=0
          ELSE
            IUTAIR=1
          END IF
 601      write(6,*) ' How many colours?'
          read(5,*,END=999,ERR=601) NCD
C         IF (NCO.EQ.0) GO TO 221
    	  call setairmass(amt,idco,color,ncd)
          write(6,150) (I,COLOR(I),I=1,NCO)
 150    FORMAT(' Colour CODES for standard data are: ',
     *             5('#',I1,1X,A1,2X))
c     *             <NCO>('#',I1,1X,A1,2X))
 221      write(6,*) ' ENTER colour codes for input data ',
     *              'in desired order:'
          read(5,*,ERR=221,END=999) (IDCO(I),I=1,NCD)
        END IF
c
C    read data
200     NSD=0
        ZTUT=0.
        RLAT=RLAT/57.2958
	istatus=encode(ncdstr,float(ncd),1,0,idum)
c 140    FORMAT(1X,A10,<NCD>(1X,F7.3,F6.2,2X),I3,I3,F5.1)
 	fstr140='(1X,A10,'//ncdstr//'(1X,F7.3,F6.2,2X),I3,I3,F5.1)'
c 143    FORMAT(1X,A10,<NCD>(1X,F7.3,F6.2,2X),F6.3)
	fstr143='(1X,A10,'//ncdstr//'(1X,F7.3,F6.2,2X),F6.3)'
        DO 100 I=1,200
          NSD=NSD+1
          IF (IPAR.EQ.0) THEN
            IF (IDAY.NE.0) THEN
              READ(28,fstr140,END=151)STN(NSD),(STMAG(NSD,J),TS(NSD,J)
     *           ,J=1,NCD),  IOH(NSD),IOM(NSD),OS(NSD)
            ELSE
              READ(28,fstr143,END=151)STN(NSD),(STMAG(NSD,J),TS(NSD,J)
     *           ,J=1,NCD), AIRMASS(NSD)
            END IF
          ELSE
            write(6,5) nsd
 5	format(' Object',i3,': enter name (END to end): ',$)
            read(5,13,END=999) STN(NSD)
 13	format(a10)
            IF (STN(NSD).EQ.'END') GO TO 129
          END IF
C  check against standard list
          DO 125 J=1,NOS
            IF (STNAM(J).EQ.STN(NSD)) GO TO 126
 125      CONTINUE
 126      IF (J.EQ.(NOS+1)) THEN   
c  star is not in list
            write(6,*) ' ****Object ',STN(NSD),' NOT in standard '
     *             ,'catalogue, NOT ENTERED'
            NSD=NSD-1
            GO TO 100
          END IF
c  star is legit
          IDST(NSD)=J
          IF (IPAR.EQ.1) THEN                  
c read from terminal
 110        write(6,*) '   Enter Inst. Mag for ',  NCD,' filtersa:'
            read(5,*,ERR=110,END=999) (STMAG(NSD,J),J=1,NCD)
 650        write(6,21) ncd
 21	format('  Enter',i2,
     *          ' integration times (sec) <cr>=same as last:')
	    read(5,'(a60)',end=999) tstr
	    if (tstr.eq.' ') then
	      if (nsd.eq.1) then
	        go to 650
	      else
	        do 170 j=1,ncd
	          ts(nsd,j)=ts(nsd-1,j)
 170	        continue
	      end if
	    else
              read(tstr,*,ERR=650,END=999) (TS(NSD,J),J=1,NCD)
	    end if
            IF (FILEN.EQ.'t') THEN
	      if (nsd.eq.1) then
	        oh=0
	        om=0
	        os(nsd)=0
	      else
	        oh=ioh(nsd-1)
	        om=iom(nsd-1)
	        os(nsd)=os(nsd-1)
	      end if
 605          write(6,6) int(oh),int(om),os(nsd)
 6	format(' Enter UT of observation [',2i3,f5.1,']: ',$)
	      call readr4(3,oh,om,os,*999,*605,*2000,*2000)
c             read(5,*,ERR=605,END=999) IOH(NSD),IOM(NSD),OS(NSD)
 2000	      ioh(nsd)=oh
	      iom(nsd)=om
            ELSE
	      if (nsd.eq.1) then
	        airmass(1)=0.
	      else
	        airmass(nsd)=airmass(nsd-1)
	      end if
 606          write(6,7) airmass(nsd)
 7	format(' Enter airmass [',f7.3,']: ',$)
              call read1r(airmass(nsd),*999,*606)
c              read(5,*,ERR=606,END=999) AIRMASS(NSD)
              GO TO 100
            END IF
          END IF
C       DO 101 J=1,NCD
C101    STMAG(I,J)=STMAG(I,J)+2.5*ALOG10(TS(NSD,J))
C
C    Calculate airmass
118       IF (IUTAIR.EQ.1) GO TO 100
          UTH=IOH(NSD)+ZTUT
          IF (UTH.GE.24.) UTH=UTH-24.
          IF (UTH.LT.0.) UTH=UTH+24.
          UTM=IOM(NSD)
c  epoch: temp fix, all SA stars from Landolt are assumed to be
c        epoch 2000, all others, 1950.
	  if (stn(nsd)(1:2).eq.'SA') then
	    epoch=100.
	  else 
	    epoch=50.
	  end if
c precess alpha
          PREALPH=ALPH(IDST(NSD))+(FLOAT(IYR)-epoch+(MON-1)*30./365.)*
     *        (3.07+1.34*SIN(ALPH(IDST(NSD))*15./206264.8)*
     *          TAN(DEC(IDST(NSD))/206264.8))
          CALL SIDTIM(UTH,UTM,OS(NSD),IDAY,MON,IYR,RLONG,STH,STM,STS)
          write(6,*) STH,STM,STS,PREALPH/3600,1900+epoch 
          HA(I)=(STH*3600.+STM*60.+STS-PREALPH)/3600.
          IF (HA(I).GT.12.) HA(I)=HA(I)-24.
          IF (HA(I).LT.-12.) HA(I)=HA(I)+24.
          HH=HA(I)*15./57.2958
c precess dec and divide to get radian
          DD=(DEC(IDST(NSD))+20*cos(hh))/206264.8
          AIRMASS(NSD)=ABS(1./(SIN(RLAT)*SIN(DD)+COS(RLAT)*
     *            COS(DD)*COS(HH)))
 100    CONTINUE
 129	IF (IPAR.EQ.0) THEN
 151      CLOSE(UNIT=28)
          NSD=NSD-1
          RETURN1
        ELSE     
c  save data on disk
          write(6,*) ' ENTER FILE NAME TO SAVE DATA FILE ON DISK',
     *           ' (n FOR NO SAVE)'
          read(5,11,END=999) FILEN
          NSD=NSD-1
          IF (FILEN.EQ.'n') RETURN1
          OPEN(UNIT=29,FILE=FILEN,status='NEW')
c 141    FORMAT(1X,I1,2X,<NCD>(I1,1X),I3,I3,I5)
	write(6,*) '  Enter comment:'
	read(5,19) comstr
	write(29,191) comstr
 191	format(1x,a60)
	fstr141='(1X,I1,2X,'//ncdstr//'(I1,1X),I3,I3,I5)'
          WRITE(29,fstr141) NCD,(IDCO(I),I=1,NCD),IDAY,MON,IYR
          IF (IUTAIR.EQ.0) THEN
            DO 139 N=1,NSD
 139          WRITE(29,fstr140) STN(N),(STMAG(N,J),TS(N,J),J=1,NCD)
     *             ,IOH(N),IOM(N),OS(N)
          ELSE
            DO 176 N=1,NSD
 176          WRITE(29,fstr143) STN(N),(STMAG(N,J),TS(N,J),J=1,NCD)
     *                  ,AIRMASS(N)
          END IF
          CLOSE(UNIT=29)
        END IF
 999    RETURN1
        END
c////////////////////////////////////////////////////
c  subroutine to set up default airmass terms
c
	subroutine  setairmass(amt,idco,color,ncd)
	dimension amt(*),idco(*)
	character*1 color(*)
	do 100 i=1,ncd
	  if (color(idco(i)).eq.'r'.or.color(idco(i)).eq.'R') then
	    amt(i)=-0.09
	  else if (color(idco(i)).eq.'i'.or.color(idco(i)).eq.'I') then
	    amt(i)=-0.04
	  else if (color(idco(i)).eq.'g'.or.color(idco(i)).eq.'V') then
	    amt(i)=-0.12
	  else if (color(idco(i)).eq.'B') then
	    amt(i)=-0.15
	  else if (color(idco(i)).eq.'z') then
	    amt(i)=-0.03
	  else if (color(idco(i)).eq.'U') then
	    amt(i)=-0.22
	  end if
 100	continue
	return
	end
C////////////////////////////////////////////////////////////////////
c  subroutine to type out standard file
c
        SUBROUTINE TSTF(*,IPAR)
        CHARACTER*12 STNAM
        CHARACTER COLOR*1,ncostr*1,fstr20*60
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *           ,DEC(700),NOS,NCO
        IDEV=6
	istatus=encode(ncostr,float(nco),1,0,idum)
        IF (IPAR.EQ.1) IDEV=8
        WRITE(IDEV,*) '  STANDARD STAR FILE'
	fstr20='(/,9X,''NAME'','//ncostr//'(4X,A1,3X))'
        WRITE(IDEV,fstr20) (COLOR(I),I=1,NCO)
c 10     FORMAT(/,9X,'NAME',<NCO>(4X,A1,3X))
c  20    FORMAT(1X,I2,1X,A10,<NCO>(2X,F6.3),5X,F9.5,2X,F9.5)
	fstr20='(1X,I2,1X,A10,'//ncostr//'(2X,F6.3),5X,F9.5,2X,F9.5)'
        DO 100 I=1,NOS
 100    WRITE(IDEV,fstr20) I,STNAM(I),(RMAG(I,J),J=1,NCO),
     *                 ALPH(I)/3600.,DEC(I)/3600.
        RETURN1
        END
C/////////////////////////////////////////////////////////////
c  subroutine to type out standard data
C
        SUBROUTINE TDAT(*,IPAR)
        CHARACTER STNAM*12,STN*12,FILEN*18,ncdstr*1,fstr*60
        CHARACTER COLOR*1,obse*4
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *           ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *         TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
	istau=encode(ncdstr,float(ncd),1,0,idum)
c 10     FORMAT(/11X,'NAME',<NCD>(5X,A1,3X,4HTIME))
	fstr='(/11X,''NAME'','//ncdstr//'(5X,A1,3X,4HTIME))'
        IDEV=6
        IF (IPAR.EQ.1) IDEV=8
        WRITE(IDEV,*) ' STANDARD STAR OBSERVATIONS'
        WRITE(IDEV,fstr) (COLOR(IDCO(I)),I=1,NCD)
c 20     FORMAT(1X,I2,1X,A10,<NCD>(2X,F7.3,F5.1),5X,F7.3,4X,F7.3)
	fstr='(1X,I2,1X,A10,'//ncdstr//'(2X,F7.3,F5.1),5X,F7.3,4X,F7.3)'
        DO 100 I=1,NSD
 100    WRITE(IDEV,fstr) I,STN(IDST(I)),(STMAG(I,J),TS(I,J),J=1,NCD)
     *                 ,HA(I),AIRMASS(I)
        RETURN1
        END
C///////////////////////////////////////////////////////////////////////////
C
C TEST PROGRAM FOR SIDTIM
C
        SUBROUTINE STIME(*)
 100    write(6,*) ' ENTER UT DAY MONTH YR'
        read(5,*) IDAY,MON,IYR
        IF (IYR.EQ.0) RETURN1
        write(6,*) ' ENTER UT HOUR MIN SEC'
        read(5,*) IUTH,IUTM,UTS
        write(6,*) 'ENTER WEST LOGITUDE IN DEGREES'
        read(5,*) RLONG
        UTH=FLOAT(IUTH)
        UTM=FLOAT(IUTM)
        write(6,*) UTH,UTM,UTS
        CALL SIDTIM(UTH,UTM,UTS,IDAY,MON,IYR,RLONG,STH,STM,STS)
        write(6,*) STH,STM,STS
        GO TO 100
        END
c////////////////////////////////////////////// 
C PROGRAM TO CALCULATE SIDEREAL TIME GIVEN UT, DATE, AND LONGITUDE
C
        SUBROUTINE SIDTIM(UTH,UTM,UTS,IDAY,MON,IYR,RLONG,
     *                   STH,STM,STS)
        REAL*8 RU,TU,UT,STT,DJUL
        INTEGER*2 MDAY(12)
        DATA MDAY/0,31,59,90,120,151,181,212,243,273,304,334/
C  CALCULATE ULIAN DAYS SINCE 1900.00 12 UT
        DJUL=365.*IYR+FLOAT((IYR-1)/4)
        DJUL=DJUL+MDAY(MON)+IDAY
        UT=UTH*3600.+UTM*60.+UTS
        TU=(DJUL+(UT-43200.)/86400.)/36525.
        DJUL=DJUL+2415019.5
        write(6,*) 'JULIAN DATE =',DJUL
        RU=67125.836+(8640184.542+.0929*TU)*TU
        STT=UT+43200.+RU-RLONG*240.
        STH=FLOAT(INT(STT/3600.))
        STM=FLOAT(INT((STT-STH*3600.)/60.))
        STS=STT-STH*3600.-STM*60.
        STH=STH-INT(STH/24.)*24.
        RETURN
        END
c////////////////////////////////////////////////////////////
C
C   SUBROUTINE FOR STANDARD REDUCTION
C       IPAR=0 : A0 TERM ONLY
C       IPAR=1 : COLOUR TERM ONLY
C       IPAR=2 : AIRMASS TERM ONLY
C       IPAR=3 : AIRMASS + COLOUR TERMS
C
c   Nov 12 1997, definition of jc1 and jc2 changed from 
c     the actual observed filter number to the standard
c     table filter number
        SUBROUTINE REDUCE(*,IPAR,IPRINT,AMT,CT)
        DIMENSION SIGMAY(200,8),M(200,8),YFIT(200,8)
     *  ,IW(200,8),A(10,8),SIGMAA(10,8),SIGMA0(8),R(10,8),
     *     RMUL(8),PMAG(200,8),ICC(8),AMT(8),CT(8)
        CHARACTER STNAM*12,STN*12,FILEN*18,NAM*10,ncdstr*1,fstr*70
        CHARACTER obse*4,COLOR*1,filnam*18
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *           ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *        TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
        COMMON/FUNN/ X1(700),X2(700)
	common/plotdata/xx(200,8),yy(200,8),nps(8),a0(8),a1(8),a2(8),
     *                  itype,numo(200,8)
	itype=ipar
        NTERMS=IPAR
        IF (IPAR.GE.2) NTERMS=IPAR-1
	istatus=encode(ncdstr,float(ncd),1,0,idum)
        DO 200 I=1,NCD
c  determine default color term
          JC1(I)=idco(I+1)
          JC2(I)=idco(I)
c input default airmass terms
c	  if (color(i).eq.'I'.or.color(i).eq.'i') then
c	     amt(i)=-.037
c	  else if (color(i).eq.'R'.or.color(i).eq.'r') then
c	     amt(i)=-.052
c	  else if (color(i).eq.'V') then
c	     amt(i)=-.10
c	  else if (color(i).eq.'g') then
c	     amt(i)=-.11
c	  else if (color(i).eq.'B') then
c	     amt(i)=-.18 
c	  else if (color(i).eq.'U') then
c	     amt(i)=-.30 
c	  end if
          DO 201 J=1,NSD
 	    IW(J,I)=1
 201	  continue
 200	continue
        JC1(NCD)=idco(NCD-1)
        CALL TDAT(*257,0)
 257     IF (IPRINT.EQ.1) THEN
	  filnam='phored.out'
	  call opennew(8,1,18,filnam,*999)
          WRITE(8,*)
          WRITE(8,*)
          WRITE(8,*) ' ////////////////////////////////////////////////'
          CALL TDAT(*258,1)
        END IF
c 55     FORMAT(1X,' DEFAULT COLOURS:',/,2X,<NCD>(A1,'-',A1,2X))
 258	fstr='(1X,'' DEFAULT COLOURS:'',/,2X,'//ncdstr//
     *	           '(A1,''-'',A1,2X))'
        write(6,fstr)(COLOR(JC1(J)),COLOR(JC2(J)),J=1,NCD)
        write(6,21) 
 21	format(' <Cr> To continue, enter any key for changes: ',$)
        read(5,78,end=999) NAM
 78     FORMAT(A10)
        IF (NAM.EQ.' ') GO TO 220
C
C   DELETE UNWANTED OBJECTS FROM LIST
C
 77     write(6,22) 
 22	format(' Do you want to delete any objects? y/n ',$)
        read(5,78,end=999) NAM
        IF (NAM.EQ.'y') THEN
 600       write(6,1) 
 1	format(' Enter obj number, 0 TO END: ',$)
	  call read1i(nn,*999,*600)
          IF (NN.GT.NSD) GO TO 600
          IF (NN.EQ.0) GO TO 179
 601      write(6,74) (J,COLOR(idco(j)),J=1,ncd)
 74     FORMAT(' Enter filter number(s): ',8('#',I1,1X,A1,', '))
	  fstr='('//ncdstr//'(I1,1X))'
          read(5,fstr,ERR=601,END=167) (ICC(J),J=1,NCD)
c  change weight to 0
 167      DO 168 L=1,J-1
	    IW(NN,ICC(L))=0
 168	  continue
          GO TO 600
        END IF
c
C change colour term
 179    write(6,2) 
 2	format(' Do you want to change reference colors? y/n: ',$)
        read(5,78) NAM
        IF (NAM.EQ.'y') THEN
          write(6,40) (COLOR(J),J,J=1,nco)
 40     FORMAT(2X,'COLOUR CODE: ',8(A1,'=',I1,', '))
          DO 146 L=1,NCD
 602         write(6,41) COLOR(IDCO(L)),jc1(l),jc2(l)
 41     FORMAT(2X,' Color term (2 integers) for ',A1,
     *     ' [',2i2,']: ',$)
	    call read2i(jc1(l),jc2(l),*999,*602)
            write(6,25) COLOR(JC1(L)),COLOR(JC2(L))
 25     FORMAT(10X,' Color term is:',2X,A1,'-',A1)
 146	  continue
        END IF
 45     IF (IPAR.EQ.3) GO TO 220
C
C  SET DEFAULT AIRMASS AND COLOUR TERMS
        IF (IPAR.LE.1) write(6,210) (AMT(J),J=1,NCD)
 210    FORMAT(5X,'DEFAULT AIRMASS COEF''S',/,10X,8(F6.3,3X))
        IF (IPAR.EQ.0.OR.IPAR.EQ.2) write(6,211) (CT(J),J=1,NCD)
 211    FORMAT(5X,'DEFAULT COLOUR COEF''S',/,10X,8(F6.3,3X))
        IF (IPAR.LE.1) THEN
          write(6,3) 
 3	format(' Do you want to change airmass? y/n [n]: ',$)
          read(5,78,END=999) NAM
          IF (NAM.EQ.'y') THEN
212         write(6,*) ' ENTER',NCD,' NEW AIRMASS COEF''S:'
            read(5,*,ERR=212,END=212) (AMT(J),J=1,NCD)
          END IF
        END IF
        IF (IPAR.EQ.0.OR.IPAR.EQ.2) THEN
          write(6,4) 
 4	format(' Do you want to change colour coefficients? y/n [n]: ',$)
          read(5,78,END=999) NAM
          IF (NAM.EQ.'y') THEN
 213        write(6,*) ' ENTER',NCD,' NEW COLOUR COEFF''S:'
            read(5,*,ERR=213) (CT(J),J=1,NCD)
          END IF
 220    END IF
C
C  SET AIRMASS AND COLOUR TERMS TO 0.0 IF THEY ARE TO BE DETERMINED
        IF (IPAR.EQ.1.OR.IPAR.EQ.3) THEN
          DO 222 J=1,8
 222        CT(J)=0.
        END IF
        IF (IPAR.EQ.2.OR.IPAR.EQ.3) THEN
          DO 223 J=1,8
 223        AMT(J)=0.
        END IF
C
C  REDUCTION FOR EACH COLOUR
c
        DO 101 J=1,NCD
          NPS(j)=0
          DO 100 I=1,NSD
c setting weight to 0 for those stars without data in both
c    filters
            IF (abs(STMAG(I,J)).lt.0.00001) IW(I,J)=0
            IF (RMAG(IDST(I),JC1(J)).EQ.0.) IW(I,J)=0
            IF (RMAG(IDST(I),JC2(J)).EQ.0.) IW(I,J)=0
            IF (IW(I,J).EQ.0) GO TO 100
            NPS(j)=NPS(j)+1
	    numo(nps(j),j)=i
C       TYPE*,RMAG(IDST(I),IDCO(J)),IDST(I),IDCO(J),IW(I,J)
C
C DEPENDENT VARIABLE X = COLOUR IF COLOUR TERM IF CT ARE TO BE DETERMINED
C                      = AIRMASS IF ONLY AMT ARE TO BE DETERMINED
            xx(NPS(j),j)=RMAG(IDST(I),JC1(J))-RMAG(IDST(I),JC2(J))
C   PMAG= INS MAG PER SECOND, CORRECTED FOR DEFAULT TERMS, ct and
c            amt have been  set to 0's if they are to be solved
            PMAG(I,J)=STMAG(I,J)+2.5*ALOG10(TS(I,J))+CT(J)*xx(NPS(j),j)
     *             +AMT(J)*AIRMASS(I)
c  x1 and x2 are common blocks giving the airmass and color 
c     in the fitting function
c         x1=color or airmass if fitting only color or airmass
c         x2=color (when fitting both airmass and color)
c         y=actual mag - instru mag corrected for default terms 
            X1(NPS(j))=AIRMASS(I)
            X2(NPS(j))=xx(NPS(j),j)
            IF (IPAR.EQ.1) X1(NPS(j))=xx(NPS(j),j)
            yy(NPS(j),j)=RMAG(IDST(I),IDCO(J))-PMAG(I,J)
 100      CONTINUE
          IF (IPAR.ne.0) then
C
C   SOLVE FOR A0, A1 AND/OR A2
C
            CALL MLRF(xx(1,j),yy(1,j),SIGMAY(1,J),NPS(j),NTERMS,M(1,J),
     *        0,YFIT(1,J),A0(J),A(1,J),SIGMA0(J),SIGMAA(1,J),R(1,J),
     *            RMUL(J),CHISQ,FTEST)
          else     
C
C  FOR BOTH COLOUR AND AIRMASS TERMS = DEFAULT
C      I.E. SOLVE FOR A0 ONLY
            A0(J)=0.
            DO 233 L=1,NPS(j)
              A0(J)=A0(J)+yy(L,j)
 233	    continue
            A0(J)=A0(J)/FLOAT(NPS(j))
            DO 234 L=1,NPS(j)
              YFIT(L,J)=A0(J)
 234	    continue
	  end if
C
C  OUT PUT RESULTS
C
 231      write(6,10) COLOR(IDCO(J)),IDAY,MON,IYR,
     *       COLOR(JC1(J)),COLOR(JC2(J))
          IF (IPAR.EQ.2) write(6,5)
          IF (IPAR.EQ.1) write(6,6)
          IF (IPAR.EQ.0) write(6,7)
          write(6,12)
 7      FORMAT(15X, '//Default airmass and color terms used//')
 6      FORMAT(15X, '//DEFAULT airmass terms used//')
 5      FORMAT(15X, '//DEFAULT color termss used//')
          IF (IPRINT.eq.1) then
            WRITE(8,*)
            WRITE(8,*)
            WRITE(8,10) COLOR(IDCO(J)),IDAY,MON,IYR
            WRITE(8,25) COLOR(JC1(J)),COLOR(JC2(J))
            IF (IPAR.EQ.2) WRITE(8,2)
            IF (IPAR.EQ.1) WRITE(8,1)
            IF (IPAR.EQ.0) WRITE(8,7)
            WRITE(8,12)
	  end if
 	  SUM=0.
          NP=0
          DO 102 I=1,NSD
            clr=RMAG(IDST(I),JC1(J))-RMAG(IDST(I),JC2(J))
            IF (IW(I,J).EQ.0) GO TO 750
            NP=NP+1
            P=RMAG(IDST(I),IDCO(J))
            Q=YFIT(NP,J)+PMAG(I,J)
            SUM=SUM+(P-Q)**2
            GO TO 751
 750        P=0.
            Q=0.
 751        IF (IPRINT.EQ.1) 
     *       write(8,11) I,STN(I),PMAG(I,J),P,Q,P-Q,clr,airmass(i)
 	    write(6,11) I,STN(I),PMAG(I,J),P,Q,P-Q,clr,airmass(i)
 102	  continue
 11     FORMAT(1X,I2,1X,A10,2X,F7.3,2X,F6.3,3X,F6.3,3X,F6.3,
     *      5x,f6.3,2x,f5.3)
 12     FORMAT('     NAME',6X,'MAG(INST)',1X,'STD MAG',
     *         2X,'CAL MAG',2X,'RESIDUAL',4x,'color',2x,'airmass')
          SUM=SQRT(SUM/(NP-1))
          A1(j)=AMT(J)
          A2(j)=CT(J)
	  istatus=encode(ncdstr,float(nterms+1),1,0,idum)
	  fstr='(18X,'//ncdstr//'(F10.3),8X,''r.m.s.:'',F9.4)'
          IF (IPAR.EQ.1) A2(j)=A(1,J)
          IF (IPAR.EQ.2.OR.IPAR.EQ.3) A1(j)=A(1,J)
          IF (IPAR.EQ.3) A2(j)=A(2,J)
          IF (IPRINT.EQ.1) then
            WRITE(8,13) A0(J),A1(j),A2(j)
            IF (IPAR.GT.0)
     *         WRITE(8,14) SIGMA0(J),(SIGMAA(I,J),I=1,NTERMS)
c 15     FORMAT(18X,<NTERMS+1>(F10.3),8X,'r.m.s.:',F9.4)
            WRITE(8,fstr) RMUL(J),(R(I,J),I=1,NTERMS),SUM
	  end if
 701      write(6,13) A0(J),A1(j),A2(j)
          IF (IPAR.GT.0)
     *         write(6,14) SIGMA0(J),(SIGMAA(I,J),I=1,NTERMS)
 	  write(6,fstr) RMUL(J),(R(I,J),I=1,NTERMS),SUM
 101	continue
	if (iprint.eq.1) close(unit=8)
 10     FORMAT(/' Calibration for filter: ',A1,' date: '
     *          ,I2,2('/',I2),3X,' Color term: ',A1,'-',A1)
 13     FORMAT(/2X,'MAG - MAG(INST) = ',F7.3,' + ',F6.3,
     *         ' *AIRMASS + ',F6.3,' *COLOR')
 14     FORMAT(17X,F7.5,3X,F7.5,3X,F7.5)
 999    RETURN1
        END
C/////////////////////////////////////////////////////////////
C
        FUNCTION FCTN(X,I,J,M)
        DIMENSION X(1),M(1)
        COMMON/FUNN/ X1(700),X2(700)
        IF (J.EQ.1) FCTN=X1(I)
        IF (J.EQ.2) FCTN=X2(I)
        RETURN
        END
c//////////////////////////////////////////////////////
c subroutine to plot fits
c
	subroutine plotdel(*)
	character splot*18,hplot*18,dev*18,xlabel*20,ylabel*20,
     *      tlabel*20, color*1,obse*4,stnam*12,stn*12
	dimension tx(2),ty(2)
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *           ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *        TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
	common/plotdata/xx(200,8),yy(200,8),nps(8),a0(8),a1(8),a2(8),
     *                  itype,numo(200,8)
        common/plotdev/hplot,splot,ipinit
	if (splot.ne.'/xwin'.and.splot.ne.'/tek') then
	  write(6,1) splot
 1	format('  Enter plotting device [',a18,']: ',$)
	  read(5,'(a18)',end=999) splot
	end if
 600	write(6,2) ncol
 2	format('  Enter filter number for plotting [',i1,']: ',$)
	call read1i(ncol,*999,*600)
	xmint=1000
	xmaxt=-1000
	ymint=1000
	ymaxt=-1000
	do 100 i=1,nps(ncol)
	  if (xx(i,ncol).lt.xmint) xmint=xx(i,ncol)
	  if (xx(i,ncol).gt.xmaxt) xmaxt=xx(i,ncol)
	  if (yy(i,ncol).lt.ymint) ymint=yy(i,ncol)
	  if (yy(i,ncol).gt.ymaxt) ymaxt=yy(i,ncol)
 100	continue
	xmint=xmint-0.1
	xmaxt=xmaxt+0.1
	ymint=ymint-0.1
	ymaxt=ymaxt+0.1
	ylabel='Delta mag'
	if (itype.eq.1.or.itype.eq.3.or.itype.eq.0) then
	    xlabel=COLOR(JC1(ncol))//'-'//COLOR(JC2(ncol))
	else
	  xlabel='air mass'
	end if
	tlabel='filter '//color(ncol)
	dev=splot
 770	call openpg(dev,ipinit,*999)
        CALL pgsch(1.)
 790    call pgsci(1)
	CALL PGENV(xmint,xmaxt,YMINT,YMAXT,0,0)
        CALL PGLABEL(XLABEL,YLABEL,TLABEL)
	call pgsci(5)
	write(6,*) ' # of pts: ',nps(ncol)
	CALL PGPOINT(nps(ncol),xx(1,ncol),yy(1,ncol),3)
	call pgsci(2)
	tx(1)=-10
	tx(2)=10
	if (itype.eq.1.or.itype.eq.3) then
	  ty(1)=-10*a2(ncol)+a0(ncol)
	  ty(2)=10*a2(ncol)+a0(ncol)
	else if (itype.eq.0) then
c for fitting only 0 point, here the delta mag yy values have
c  already the color and airmass terms taken out, so the effective
c   slope should be 0
	  ty(1)=a0(ncol)
	  ty(2)=a0(ncol)
	else
	  ty(1)=-10*a1(ncol)+a0(ncol)
	  ty(2)=10*a1(ncol)+a0(ncol)
	end if
	call pgline(2,tx,ty)
 799	CALL ITACUR(xmint,xmaxt,YMINT,YMAXT,CX,CY,DEV,ipinit,IRET)
        IF (IRET.EQ.1) THEN
            GO TO 790
        ELSE IF (IRET.EQ.2) THEN
            if (dev.ne.'/sun'.and.dev.ne.'/xwin') CALL closepg(ipinit)
            GO TO 770
        ELSE IF (IRET.EQ.-1) THEN
            write(idout,*) CX,CY
            GO TO 799
        ELSE IF (IRET.EQ.0) THEN
            if (dev.ne.'/sun'.and.dev.ne.'/xwin') CALL closepg(ipinit)
            RETURN1
        ELSE IF (IRET.EQ.3.OR.IRET.EQ.4) THEN
            GO TO 790
        ELSE IF (IRET.EQ.10) THEN
            GO TO 790
        END IF
 999	return1
	end
C/////////////////////////////////////////////////////////////
C
C   SUBROUTINE FOR MULTIPLE LINEAR REGRESSION
C
        SUBROUTINE MLRF(X,Y,SIGMAY,NPTS,NTERMS,M,MODE,YFIT,
     #                 A0,A,SIGMA0,SIGMAA,R,RMUL,CHISQR,FTEST)
        REAL*8 ARRAY,SUM,YMEAN,SIGMA,CHISQ,XMEAN,SIGMAX
        DIMENSION X(1),Y(1),SIGMAY(1),M(1),YFIT(1),A(1),SIGMAA(1)
     #           , R(1)
        DIMENSION WEIGHT(700),XMEAN(700),SIGMAX(10),ARRAY(10,10)
        COMMON/FUNN/ X1(700),X2(700)
C
 11     SUM=0.
        YMEAN=0.
        SIGMA=0.
        CHISQ=0.
        RMUL=0.
        DO 17 I=1,NPTS
 17     YFIT(I)=0.
 21     DO 28 J=1,NTERMS
        XMEAN(J)=0.
        SIGMAX(J)=0.
        R(J)=0.
        A(J)=0.
        SIGMAA(J)=0.
        DO 28 K=1,NTERMS
 28     ARRAY(J,K)=0.
C
 30     DO 50 I=1,NPTS
 31     IF (MODE) 32,37,39
 32     IF (Y(I)) 35,37,33
 33     WEIGHT(I)=1./Y(I)
        GO TO 41
 35     WEIGHT(I)=1./(-Y(I))
        GO TO 41
 37     WEIGHT(I)=1.
        GO TO 41
 39     WEIGHT(I)=1./SIGMAY(I)**2.
 41     SUM=SUM+WEIGHT(I)
        YMEAN=YMEAN+WEIGHT(I)*Y(I)
        DO 44   J=1,NTERMS
 44     XMEAN(J)=XMEAN(J)+WEIGHT(I)*FCTN(X,I,J,M)
 50     CONTINUE
 51     YMEAN=YMEAN/SUM
        DO 53 J=1,NTERMS
 53     XMEAN(J)=XMEAN(J)/SUM
        FNPTS=NPTS
        WMEAN=SUM/FNPTS
        DO 57 I=1,NPTS
 57     WEIGHT(I)=WEIGHT(I)/WMEAN
C
 61     DO 67 I=1,NPTS
        SIGMA=SIGMA+WEIGHT(I)*(Y(I)-YMEAN)**2.
        DO 67 J=1,NTERMS
        XF=FCTN(X,I,J,M)-XMEAN(J)
        SIGMAX(J)=SIGMAX(J)+WEIGHT(I)*(XF)**2
        R(J)=R(J)+WEIGHT(I)*(XF)*(Y(I)-YMEAN)
        DO 67 K=1,J
 67     ARRAY(J,K)=ARRAY(J,K)+WEIGHT(I)*(XF)*(FCTN(X,I,K,M)-XMEAN(K))
 71     FREE1=NPTS-1
 72     SIGMA=DSQRT(SIGMA/FREE1)
        DO 78 J=1,NTERMS
 74     SIGMAX(J)=DSQRT(SIGMAX(J)/FREE1)
        R(J)=R(J)/(FREE1*SIGMAX(J)*SIGMA)
        DO 78 K=1,J
        ARRAY(J,K)=ARRAY(J,K)/(FREE1*SIGMAX(J)*SIGMAX(K))
 78     ARRAY(K,J)=ARRAY(J,K)
C
 81     CALL MATINV(ARRAY,NTERMS,DET)
        IF (DET) 101,91,101
 91     A0=0.
        SIGMA0=0.
        RMUL=0.
        CHISQ=0.
        FTEST=0.
        GO TO 150
C
 101    A0=YMEAN
 102    DO 108 J=1,NTERMS
        DO 104 K=1,NTERMS
 104    A(J)=A(J)+R(K)*ARRAY(J,K)
 105    A(J)=A(J)*SIGMA/SIGMAX(J)
 106    A0=A0-A(J)*XMEAN(J)
 107    DO 108 I=1,NPTS
 108    YFIT(I)=YFIT(I)+A(J)*FCTN(X,I,J,M)
 111    DO 113 I=1,NPTS
        YFIT(I)=YFIT(I)+A0
 113    CHISQ=CHISQ+WEIGHT(I)*(Y(I)-YFIT(I))**2
        FREEN=MAX0(1,NPTS-NTERMS-1)
 115    CHISQ=CHISQ*WMEAN/FREEN
C
 121    IF (MODE) 122, 124, 122
 122    VARNCE=1./WMEAN
        GOTO 131
 124    VARNCE=CHISQ
 131    DO 133 J=1,NTERMS
 132    SIGMAA(J)=ARRAY(J,J)*VARNCE/(FREE1*SIGMAX(J)**2)
        SIGMAA(J)=SQRT(SIGMAA(J))
 133    RMUL=RMUL+A(J)*R(J)*SIGMAX(J)/SIGMA
        FREEJ=NTERMS
 135    FTEST=(RMUL/FREEJ)/((1.-RMUL)/FREEN)
        IF (RMUL.LT.0.) RMUL=0.
 136    RMUL=SQRT(RMUL)
 141    SIGMA0=VARNCE/FNPTS
        DO 145 J=1,NTERMS
        DO 145 K=1,NTERMS
 145    SIGMA0=SIGMA0+VARNCE*XMEAN(J)*XMEAN(K)*ARRAY(J,K)/
     *          (FREE1*SIGMAX(J)*SIGMAX(K))
 146    SIGMA0=SQRT(SIGMA0)
 150    RETURN
        END
C
C
        SUBROUTINE MATINV(ARRAY,NORDER,DET)
        REAL*8 ARRAY,AMAX,SAVE
        DIMENSION ARRAY(10,10),IK(10),JK(10)
 10     DET=1.
 11     DO 100 K=1,NORDER
C
        AMAX=0.
 21     DO 30 I=K,NORDER
        DO 30 J=K,NORDER
 23     IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30
 24     AMAX=ARRAY(I,J)
        IK(K)=I
        JK(K)=J
 30     CONTINUE
C
 31     IF (AMAX) 41,32,41
 32     DET=0.
        GO TO 140
 41     I=IK(K)
        IF (I-K) 21,51,43
 43     DO 50 J=1,NORDER
        SAVE=ARRAY(K,J)
        ARRAY(K,J)=ARRAY(I,J)
 50     ARRAY(I,J)=-SAVE
 51     J=JK(K)
        IF (J-K) 21,61,53
 53     DO 60 I=1, NORDER
        SAVE=ARRAY(I,K)
        ARRAY(I,K)=ARRAY(I,J)
 60     ARRAY(I,J)=-SAVE
C
 61     DO 70 I=1,NORDER
        IF (I-K) 63,70,63
 63     ARRAY(I,K)=-ARRAY(I,K)/AMAX
 70     CONTINUE
 71     DO 80 I=1,NORDER
        DO 80 J=1,NORDER
        IF (I-K) 74,80,74
 74     IF (J-K) 75,80,75
 75     ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
 80     CONTINUE
 81     DO 90 J=1, NORDER
        IF (J-K) 83,90,83
 83     ARRAY(K,J)=ARRAY(K,J)/AMAX
 90     CONTINUE
        ARRAY(K,K)=1./AMAX
 100    DET=DET*AMAX
 101    DO 130 L=1,NORDER
        K=NORDER-L+1
        J=IK(K)
        IF (J-K) 111,111,105
 105    DO 110 I=1,NORDER
        SAVE=ARRAY(I,K)
        ARRAY(I,K)=-ARRAY(I,J)
 110    ARRAY(I,J)=SAVE
 111    I=JK(K)
        IF (I-K) 130,130,113
 113    DO 120 J=1,NORDER
        SAVE=ARRAY(K,J)
        ARRAY(K,J)=-ARRAY(I,J)
 120    ARRAY(I,J)=SAVE
 130    CONTINUE
 140    RETURN
        END
c////////////////////////////////////////////////////////////////
c  subroutine to read in two .cat file and do boot strap calibration
c    the data from the calibrating frame is put into the standard
c    star array, and the main data into the data array
c
c    1. It first read in two files of ppp numbers, giving option of
c     rearranging them so that they match
c    2. It read in the mag for these stars from the .cat file
c    
	subroutine calcat(*)
	parameter(nsmax=200)
	dimension npboot(nsmax),ixboot(nsmax),iyboot(nsmax),ipkboot(nsmax)
     *     ,npdat(nsmax),ixdat(nsmax),iydat(nsmax),ipkdat(nsmax)
	character pnumbootnam*24,pnumdatnam*24,bootcatnam*24,datcatnam*24,
     *      icom*1,dumstr*12
	character stnam*12,stn*12,obse*4,color*1
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *                  ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *       TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
 601	write(6,1)
 1	format(' Enter ppp number file names for boot and data images: ',$)
	read(5,*,end=999) pnumbootnam,pnumdatnam
 602	write(6,2)
 2	format(' Enter .cat names for the boot and data images',
     *      '(.cat appended): ',$)
	read(5,*,end=999) bootcatnam,datcatnam
c read file
	open(25,file=pnumbootnam,status='old',err=500)
	write(6,*) '  reading boot strap objects ',pnumbootnam
	call readpnum(25,numboot,npboot,ixboot,iyboot,ipkboot)
	close(unit=25)
	goto 501
 500	write(6,*) '  *** boot frame ppp number file does not exist **'
	goto 601
 501	open(25,file=pnumdatnam,status='old',err=502)
	write(6,*) '  reading boot strap objects ',pnumdatnam
	call readpnum(25,numdat,npdat,ixdat,iydat,ipkdat)
	close (unit=25)
	goto 503
 502	write(6,*) '  *** data frame ppp number file does not exist **'
	goto 601
c allow the user to rearrange the ppp number so that they match
 503	if (numdat.ne.numboot) then
	  write(6,*) ' # of data point ',numdat,' not eq # of cal std'
     *	           ,numboot,'  ** aborted ** '
	  return1
	end if
	call arrpppnum(numboot,npboot,ixboot,iyboot,ipkboot,'boot',*999)
  	call arrpppnum(numdat,npdat,ixdat,iydat,ipkdat,'data',*999)
c read mag from the .cat files
	call rdcat(bootcatnam,numboot,npboot,0)
	call rdcat(datcatnam,numdat,npdat,1)
c allow the user to change id of objects, or delete unwanted ones
 603	write(6,19)
 19	format('  N   bootpn datapn        I           R           B  ',
     *          '      V ')
	do 150 i=1,nos
	  write(6,20) i,npboot(i),npdat(i),rmag(i,1),stmag(i,1),
     *        rmag(i,2),stmag(i,2),rmag(i,3),stmag(i,3),
     *        rmag(i,4),stmag(i,4),ixboot(i),iyboot(i)
 20	format(i3,1x,2i6,1x,4(2f7.2,1x),1x,2i4)
 150	continue
 604	write(6,21) 
 21	format(' Enter edit command: s=switch(data), d=delete(both),t,q: '
     *         ,$)
	read(5,'(a1)',end=999) icom
	if (icom.eq.'t') then
	  goto 603
	else if (icom.eq.'q') then
	  goto 999
	else if (icom.eq.'s') then
 605	  write(6,22) is1,is2
 22	format('  Enter (seq) numbers of objs to be switched (data)[',2i3,
     *         ']: ',$)
	  call read2i(is1,is2,*999,*605)
	  call switchn(npdat(is1),npdat(is2))
	  call switchn(ixdat(is1),ixdat(is2))
	  call switchn(iydat(is1),iydat(is2))
	  call switchn(ipkdat(is1),ipkdat(is2))
	  dumstr=stn(is1)
	  stn(is1)=stn(is2)
	  stn(is2)=dumstr
	  do 140 kk=1,nco
	    call switchr(stmag(is1,kk),stmag(is2,kk))
 140	  continue
	  goto 604
	else if (icom.eq.'d') then
 606	  write(6,23) is1	
 23	format('  Enter (seq) number of obj to be deleted [',i3,']: ',$)
	  call read1i(is1,*999,*605)
	  do 141 k=is1,nos-1
	    stnam(k)=stnam(k+1)
	    do 142,i=1,nco
	      rmag(k,i)=rmag(k,i+1)
	      stmag(k,i)=stmag(k,i+1)
 142	    continue
	    stn(k)=stn(k+1)
	    idst(k)=idst(k+1)
 141	  continue
	  nos=nos-1
	  nsd=nsd-1  
	  goto 604
	end if
 999	return1
	end    
c/////////////////////////////////////////////////////
c
c  subroutine to read pppnumber file
c
	subroutine readpnum(idev,i,np,ix,iy,ipk)
	dimension np(*),ix(*),iy(*),ipk(*)
	character*8 dumstr1,dumstr2,dumstr3
c  comment line
	read(idev,*)
	do 100 i=1,100
c	  read(idev,*,end=101) np(i),dumstr1,ix(i),iy(i),dumstr2,
c    *       dumstr3,ipk(i)
	  read(idev,1,end=101) np(i),ix(i),iy(i),ipk(i)
 1	format(i5,8x,2i6,14x,i8)
 100	continue
 101	i=i-1
	write(6,*) '  Number of objects read: ',i, ' (p#,x,y,pk)'
c	do 110 k=1,i
c	  write(6,2) np(k),ix(k),iy(k),ipk(k)
c 2	format(i6,2x,2i7,3x,i7)
c110	continue
	return
	end
c//////////////////////////////////////////////////////
c
c  subroutine to rearrange ppp number array
c
	subroutine arrpppnum(num,np,ix,iy,ipk,typestr,*)
	dimension np(*),ix(*),iy(*),ipk(*)
	dimension itemp(200),indx(200)
	character typestr*4,astr*2
 601	write(6,1) typestr,num
 1	format('  Rearranging orders of ',a4,' ppp#s: ',i4/
     *    '   Enter ordering desired (pn: ppp#; ',
     *    'x/-x: x(small/large), y/-y, pk: peak): ',$)
	read(5,'(a2)',end=999) astr
	if (astr.eq.'pn') then
c  arrange in ascending order of ppp numbers
	  call nindexx(num,np,indx)
	else if (astr.eq.'x') then
	  call nindexx(num,ix,indx)
	else if (astr.eq.'y') then
	  call nindexx(num,iy,indx)
	else if (astr.eq.'ipk') then
	  call nindexx(num,ipk,indx)
	else if (astr.eq.'-x') then
	  do 100 i=1,num
	    itemp(i)=10000-ix(i)
 100	  continue
	  call nindexx(num,itemp,indx)
	else if (astr.eq.'-y') then
	  do 101 i=1,num
	    itemp(i)=10000-iy(i)
 101	  continue
	  call nindexx(num,itemp,indx)
	end if
	call rearran(num,indx,np,ix,iy,ipk)
	do 110 k=1,num
	  write(6,2) np(k),ix(k),iy(k),ipk(k)
 2	format(i6,2x,2i7,3x,i7)
 110	continue
 	return
 999	return1
	end
c///////////////////////////////////////////
c  subroutine to re-arrange the 3 arrays
c
	subroutine rearran(num,indx,np,ix,iy,ipk)
	dimension np(*),ix(*),iy(*),ipk(*),indx(*)
	dimension npt(num),ixt(num),iyt(num),ipkt(num)
	do 100 i=1,num
	  npt(i)=np(indx(i))
	  ixt(i)=ix(indx(i))
	  iyt(i)=iy(indx(i))
	  ipkt(i)=ipk(indx(i))
 100	continue
	do 101 i=1,num
	  np(i)=npt(i)
	  ix(i)=ixt(i)
	  iy(i)=iyt(i)
	  ipk(i)=ipkt(i)
 101	continue
	return
	end
c///////////////////////////////////////////////////
c  subroutine to read .cat file
c    if ipar=0 read the data into standard file
c            1                    measured star file
c
	subroutine rdcat(filnam,nump,npppnum,ipar)
	character filnam*24,filt*24,stnam*12,stn*12,obse*4,color*1
	dimension npppnum(nump)
        COMMON/STSTAR/ STNAM(700),RMAG(700,8),COLOR(8),ALPH(700)
     *                  ,DEC(700),NOS,NCO
        COMMON/STMAG/ STN(200),STMAG(200,8),IDST(200),IDCO(8),HA(200),
     *       TS(200,8),AIRMASS(200),NSD,NCD,IDAY,MON,IYR,obse,
     *        jc1(8),jc2(8)
	filt=filnam(1:lnblnk(filnam))//'.cat'
	open(unit=29,file=filt,status='old',err=999)
	if (ipar.eq.0) then
c  read standard comparison stars
	  read(29,*)
	  read(29,1) nco,nobtem
 1	format(7x,i2,8x,i5)
	  read(29,*)
c  read color info
	  do 100 i=1,nco
	    read(29,'(1x,a1)') color(i)
 100	  continue
	write(6,*) ' boot ',nco,nobtem
	  call rdcatobj(29,nump,npppnum,nco,nobtem,700,rmag)
	  nos=nump
	  do 110 i=1,nos
	    istat=encode2(stnam(i),float(npppnum(i)),1,0,idum)
	    stn(i)=stnam(i)
c	write(6,*) npppnum(i),'  ',stnam(i),stn(i)
 110	  continue
	else
c read data frame stars
	  read(29,*)
	  read(29,1) ndum,nobtem
	write(6,*) ' data ',ndum,nobtem
	  call rdcatobj(29,nump,npppnum,nco,nobtem,200,stmag)
	  nsd=nump
	  ncd=nco
	  do 101 i=1,nco
	    read(29,*) 
	    idco(i)=i
	    do 111 k=1,nsd
	      ts(k,i)=1.
 111	    continue
 101	  continue
	  do 102 i=1,nsd
	    idst(i)=i
	    airmass(i)=1.0
 102	  continue
	  idat=1
	  obse='ctio'
	end if
	close(unit=29)
	return
 999	write(6,*) ' *** no such file: ',filt
	close(unit=29)
	return
	end
c////////////////////////////////////////////////
c   reading the .cat object entries
c
	subroutine rdcatobj(idev,np,npppn,nco,ntot,ns,smag)
	dimension smag(ns,8),npppn(*)
	character str*100
	do 100 i=1,np  
	  rewind(unit=idev)
c	write(6,*) '  finding star',i,npppn(i)
	  do 101 k=1,10
	    read(29,*)
 101	  continue
	  do 120 m=1,ntot
	    read(29,'(i5)') itemp
	    if (itemp.eq.npppn(i)) then
	      backspace 29
	      read(29,'(a100)') str
	      do 102 j=1,nco
	        m1=38+(j-1)*17
	        read(str(m1:m1+5),'(f6.2)') smag(i,j)
 102	      continue
	      goto 100
	    end if
 120	  continue
	  write(6,*) '   *** warning: star not found!! ** '
 100	continue
	return
	end
c///////////////////////////////////////////////////////////
c  subroutine to switch 2 integers
c
	subroutine switchn(n1,n2)
	ndum=n1
	n1=n2
	n2=ndum
	return
	end
c///////////////////////////////////////////////////////////
c  subroutine to switch 2 real numbers
c
	subroutine switchr(r1,r2)
	rdum=r1
	r1=r2
	r2=rdum
	return
	end
c//////////////////////////////////////////////////////////
c  subroutine to convert stetson pho files to phored standard file
c
	subroutine pho2stf(*)
	character*16 phoname,outname,stname*8,stout*8
	dimension rmag(4)
 601	write(6,1)
 1	format('  Enter .pho file name (.pho appended) :',$)
	read(5,'(a16)',end=999) phoname
	outname=phoname(1:lnblnk(phoname))//'.stf'
	phoname=phoname(1:lnblnk(phoname))//'.pho'
	open(unit=25,file=phoname,status='old',err=601)
	call opennew(27,1,16,outname,*999)
c read dummy file
	read(25,*)
	do 100 i=1,1000
	  read(25,2,end=110) stname,(rmag(k),k=1,4)
 2	format(1x,a8,5x,f7.3,3(18x,f7.3))
	  if (stname(4:4).eq.'-') then
	    stout='SA'//stname(2:7)
	  else
	    stout='SA'//stname(2:4)//stname(6:8)
	  end if
	  write(27,3) stout,(rmag(k),k=4,1,-1)
 3	format(1x,'''',a8,'''',4(f7.3),' 00.000 00 00 00    00 00 00.')
 100	continue
 110	close(unit=27)
 111	close(unit=25)
 999	return1
	end