CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  STDLC_VIW.F
C    TAKES AN UNDERSAMPLED RR LYRAE LIGHT CURVE AND FITS IT WITH
C    A SET OF STANDARD LIGHT CURVES (SOLVE FOR DELTA_PHASE,
C    AMPLITUDE, INTENSITY MEAN, AND WHICH CURVE GIVES SMALLEST
C    CHI-SQUARED).  PLOT EM ALL (SM) WITH BEST CURVE BOLD.
C    THIS VERSION USES SIMPLEX OPTIMIZATION TO MINIMIZE CHI-SQUARED.
C                   ANDY LAYDEN   JAN96
C
C  R = ROBUST... SET IT SO WHEN AMOEBA DIVERGES, IT BAILS WITH CURRENT
C     BEST PARAMS, BUT SETS XSQ REAL BIG SO THAT FIT DOESN'T GET SELECTED.
C
C  THIS VERSION SPECIALLY TWEAKED FOR IHALO USE: 
C    1) POINTS PLOTTED ONCE, FITTED CURVE PLOTTED OVER -0.5 < PHASE < 1.5
C    2) EACH FITTED LC MADE IN PLACE, OBSD POINTS SHIFTED TO MATCH
C     NB: ALL INITIAL PHASES CORRECTED FOR THIS SMALL DELTA(PHI)!
C    3) OUTPUT INCLUDES FINAL FITTED AND OBSERVED (PHI-SHIFTED) DATA FILES,
C     AND 2 DATA SUMMARY FILES WHICH INCLUDE THE <V> MAGS FROM OBSD AND FIT.
C    4) THIS VERSION TWEAKED UP (WRT STDLC2.F) FOR FINAL IHALO USE!!!
C
C  THIS IS A SOFT VERSION FOR USE ON STARS WHERE AMOEBA CONVERGENCE FAILED
C   USING STDLC_IH.F.  I'VE RELAXED THE CONVERGENCE TOLERANCE AND UPPED THE
C   NUMBER OF POSSIBLE ITTNS.
C
C  THIS VERSION COPIED FROM /usr/Kracken2/layden/IHALO/PHOT/LCS/stdlc_ihr.e
C   USES ROBUST TREATMENT OF CHISQ.
C   ADAPTED FOR USE ON NGC6316 DATA, SPECIFICALLY, SELECT TEMPLATE SET.
C   ALSO INCLUDES IMAGS AND ERRORS IN OUTPUT FILE.  98dec11
C
C  THIS VERSION COPY/ADAPTED FROM 
C   /usr/Kracken3/ladyen/N6316/PERIOD/LCS/stdlc_ihr.f
C   IMPROVEMENTS INCLUDE:  DOUBLE-CHECKING ROBUSTIFICATION OF SMALL CHISQ
C   FROM BAD FITS, INCLUDING ERROR-WEIGHTING IN CHISQ TO MATCH 
C   /usr/Kracken5/layden/N6316/PERIODS/stdlc_per6a.f      ACL 98dec14
C
C  THIS VERSION COPY/ADAPTED FROM /usr/Kracken4/layden/N6316/LCS/stdlc_viw.f
C     THIS DOES THE FINAL, MONO-PERIOD FITS USING **WEIGHTED** LEAST SQRS
C     (THUS DIFFERENING FROM stdlc_finfit(2).f), AND COMPUTES <I>, 
C     V-Imin, err, Nmin FROM THE **REAL** DATA (THUS DIFFERS FROM 
C     stdlc_viw.f, AND LIKE stdlc_finfit2.f).            ACL 99jan14
C
C  THIS VERSION FROM /usr/Kracken4/layden/M54/RRL/BEN/PERIOD/LC/stdlc_finfit3.f
C     THIS DOES ALL THAT, AND CREATES A SM FILE FOR FINAL VI LC, zdat.LNN
C     ALSO ASSUMES USING stdlc.all10 TEMPLATE FILE (6RRAB, 1RRC, 1SIN,
C     1 WUMA, 1 ALGOL).                                  ACL 99jan21
C
C  THIS VERSION LIKE stdlc_finfit4.f BUT SOME BULLET-PROOFING FOR
C     NOBS=0,1 IN VIMIN ERROR COMPUTATION    ACL 99MAR01
C
C     Andrew Layden -- 1999 March 01 
C
C  THIS VERSION /usr/Kracken4/layden/M54/RRL/BEN/PERIOD/stdlc_finfit5.f
C     TIDIED UP FOR EXPORT.  2000 May 11
C
C  THIS VERSION COPIED FROM 
C     /usr/Kracken4/layden/M54/RRL/BEN/PERIOD/stdlc_export.f
C     AND ADAPTED SO INPUT MATCHES INPUT FOR stdlc_per10_export.f
C     AND THE THIS PROGRAM CAN BE USED TO CHECK APPEARANCE OF LC
C     AT CANDIDATE PERIODS AND OBTAIN FINAL LC PARAMETERS.  2000 May 24. 
C
C     compile with -e option (132 character lines): f77 stdlc_export.f -e
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       character infile*50, slcfile*50, outfile*50, datfile*50, origf*50
C  GOOD FOR UP TO 256 POINTS IN OBSD AND STD LCS, UP TO 10 STDLCS
       real v(256),p(256),ve(256),sv(10,256),sp(256)
       real imag(256),ierr(256),iimo
       real*8 hjd(256),hjdi,period,emaxo,emaxf
C  THESE ARE THE FINAL FITTED LIGHT CURVES
       real spf(256),svf(256)
C  THESE ARE THE STDLCS OF THE 4 BEST FITS, NEED ARRAYS FOR CONVENIENT STORAGE
       real slcx(4,512),slcy(4,512)
C  THESE ARE THE FINAL COEFFICIENTS, A,P0,VI; ALSO CHISQ
       real xx(3),xx1(10),xx2(10),xx3(10),xxsq(10)
       integer nptp(256),ith(10)
C       common/lcdat/nobs,p,v,nslc,sp,sv,ilc,ve
       common/lcdat1/nobs,p,v,ve,imag,ierr
       common/lcdat2/nslc,ilc,sp,sv

C   ERROR FLOOR (TO KEEP WEIGHTS REASONABLE)
	parameter (emin = 0.005)
C   OFFSET IN IMAG FOR FINAL VI LC PLOT
        parameter (offi = -0.0)
C   ERRORBAR FILTERS
        parameter (filtmag1 = 50.0, filtmag2 = -50.0, filterr = 30.0)

       write(6,90)
 90    format('Reading STD LCs from stdlc.all10')
       slcfile = 'stdlc.all10'
C       write(6,90)
C 90    format('Enter name of file w/ STD LCs (stdlc.3rr1sin2ecl or stdlc.6rrl_tmplts):  ',$)
C       read(5,100) slcfile
       write(6,91)
 91    format('Enter root name for output data files (rd.* & fd.*):  ',$)
       read(5,100) datfile
 1     write(6,80)
 80    format(//,'Enter name of input phot file (HJD,V; q=quit):  ',$)
       read(5,100) infile
       if(infile.eq.'q'.or.infile.eq.'Q') goto 2
       origf = infile
C       write(6,95)
C 95    format('Enter root name for output files:  ',$)
C       read(5,100) outfile
       outfile = infile

       write(6,96) outfile
 96    format(/,'Best fitted curve in     a',a50)
       write(6,97) outfile
 97    format('Best 4 fitted LCs in     b',a50)
       write(6,98) outfile
 98    format('Fitted coefficients in   c',a50)
       write(6,988) outfile
 988   format('4xPhi-shifted obs dat in d',a50)
       write(6,99) outfile
 99    format('4xVfit SM driver in      s',a50)
       write(6,94) outfile
 94    format('Final fitted curve in    f',a50)
       write(6,93) outfile
 93    format('Final photometry in      g',a50,/)

100    format(a)

       open(8,file=infile,status='old')
       open(9,file=slcfile,status='old')
       open(10,file='c'//outfile,status='new')
       open(11,file='a'//outfile,status='new')

C  READ IN ERROR RATIO (IE, MULTIPLY BAD-SEEING ERRORBARS BY THIS
C   VALUE, SO THEY ARE WEIGHTED **LOWER** IN THE CHISQ FIT)
	write(6,923)
923	format('Enter error ratio for bad seeing pts (eg 2.0):  ',$)
	read(5,*) erat

C  READ IN OBSERVED DATA (PHI=0 FOR POINT WITH VMAX) AND COMPUTE 
C   LC PARAMETERS DIRECTLY FROM DATA.

C  READ PERIOD FROM FIRST RECORD OF OBSERVED RRLYR FILE
C       read(8,*) period
        write(6,723) 
 723    format(/,' Enter Period of variable star:  ',$)
        read(5,*) period
C  READ DATA FROM OBSERVED RRLYR FILE (PHASE,VMAG), NOBS
       nobs = 0
       vmino = -99.
       vmaxo = 99.
       hjdmax = -9999.
       hjd0 = -9999.
       do 10 i=1,512
C          read(8,*,end=11) pi,vi,nptpi,hjdi,verri,omag,oerr
C     THIS FORMAT IS THE SAME AS USED BY stdlc_per10_export.f
          read(8,*,end=11) hjdi, nptpi, vi, verri, omag, oerr
C     CHECK FOR OUT OF BOUNDS ERROR (INDICATES NO DATA IN 1st MAG COLUMN)
          if(verri.gt.9.999) goto 10
C     IN THIS FORMAT, WE NEED TO COMPUTE A PROVISIONAL PHASE, PI,
C       BASED ON THE FIRST-ENCOUNTERED OBSERVATION (==> PI=0).
          nobs = nobs+1
          if(i.eq.1) hjd0 = hjdi
          cyc = (hjdi-hjd0)/period
          p(nobs) = cyc - float( int(cyc) )
          if(p(nobs).lt.0.0) p(nobs) = p(nobs)+1.0
          if(p(nobs).ge.1.0) p(nobs) = p(nobs)-1.0
          WRITE(6,*) NOBS, HJDI, P(NOBS), VI
C          if(pi.ge.1..or.pi.lt.0.) goto 10
C          p(nobs) = pi
          v(nobs) = vi
          nptp(nobs) = nptpi
          hjd(nobs) = hjdi
          imag(nobs) = omag
          ierr(nobs) = oerr
C   THIS PUTS A FLOOR ON ERROR, SO WEIGHTING DOESN'T GO NUTS
	  if(verri.lt.emin) verri=emin
C   THIS BIT MULTIPLIES ERROR BY ERAT FOR BAD SEEING POINTS
C    (IE, POINTS WITH NPTPI = 30,40,50...)
	  open = 10.0*(float(nptpi)/10.0 - float(nptpi/10))
	  if(abs(open).lt.0.01) verri = verri*erat
          ve(nobs) = verri
C   NOTE: VMIN/MAX IN MAG/LC SENSE (SMALL NUMBERS IS MAX).
          if(v(nobs).gt.vmino) then
            vmino = v(nobs)
            jmin = nobs
          endif
          if(v(nobs).lt.vmaxo) then
            vmaxo = v(nobs)
            jmax = nobs
          endif
          if(hjd(nobs).gt.hjdmax) lasthjd = nobs
 10       continue
 11    write(6,111) nobs
 111   format('Npts,obsd = ',i5)
       close(8)
C  COMPUTE OBSERVED AMPLITUDE
       amplo = vmino-vmaxo
C  COMPUTE RISE TIME
       WRITE(6,*) jmax,jmin,p(jmax),p(jmin)
       dphiro = p(jmax)-p(jmin)
       if(dphiro.lt.0.) dphiro = dphiro+1.
       if(dphiro.ge.1.) dphiro = dphiro-1.
C  GET INTENSITY MEAN MAG OF RAW DATA PTS
       call getimean(nobs,p,v,vimo)
       WRITE(6,*) vimo
       call getimean(nobs,p,imag,iimo)
       WRITE(6,*) iimo
C  COMPUTE EPOCH OF MAXIMUM LIGHT FROM PHASE OF LAST-OBSED DATA PT
       emaxo = hjd(lasthjd) - p(lasthjd)*period
C I MOVED THE FOLLOWING TO AFTER DATAPTS HAVE BEEN SHIFTED IN PHASE AFTER FIT
CC  COMPUTE COLOR AT MINIMUM LIGHT AND STATS.  OMIT BAD DATA PTS (EG NO I OBS)
C       cmin = 0.
C       csq = 0.
C       nmin = 0
C       do 1432 k=1,nobs
C          if(p(k).lt.0.5.or.p(k).gt.0.8) goto 1432
C          cdat = v(k)-imag(k)
C          if(cdat.lt.-5.0.or.cdat.gt.10.0) goto 1432
C          cmin = cmin + cdat
C          csq = csq + cdat*cdat
C          nmin = nmin+1
C 1432     continue
C       if(nmin.lt.1) then
C         cmino = 99.999
C         emino = 88.888
C         goto 1435
C       endif
C       cmino = cmin/float(nmin)
C       if(nmin.le.1) then
C         emino = 99.999
C         goto 1435
C       endif
C       sd = sqrt( (csq - cmin*cmin/float(nmin)) / (float(nmin-1)) )
C       emino = sd/sqrt(float(nmin))
CC  WRITE THE DERIVED DATA TO THE DATA OUTPUT FILE
C 1435  open(17,file='rd.'//datfile)
C       do 147 i=1,32000
C          read(17,100,end=148) infile
C 147       continue
C 148   write(6,4555)
C 4555  format(/,'Light Curve Parameters from Observed Data Points:')
C       write(6,455)
C 455   format('outfile    <V>o   Vmaxo  Vmino   Amplo  dPro      Emaxo     Nobs  <I>o  V-Imin  sem   Nmin')
C       write(6,456)  outfile,vimo,vmaxo,vmino,amplo,dphiro,emaxo,nobs,iimo,cmino,emino,nmin
C       write(17,456) outfile,vimo,vmaxo,vmino,amplo,dphiro,emaxo,nobs,iimo,cmino,emino,nmin
C 456   format(a10,5f7.3,f14.5,i5,3f7.3,i4)
C       close(17)


C  NOW PREPARE TO DO THE STANDARD LIGHT CURVE FITTING

C   READ STD LCS FROM FILE (phi v1 v2 v3 ... vNstd)
       read(9,*) nstd
       nslc = 0
       do 30 j=1,256
          read(9,*,end=31) sp(j),(sv(i,j),i=1,nstd)
C          if(j.le.6.or.j.ge.45) WRITE(6,*) j,sp(j),(sv(i,j),i=1,nstd)
          nslc = nslc+1
 30       continue
 31    write(6,311) nstd,nslc
 311   format('Nstd, Npts,slc = ',i2,i5,/)
       close(9)

C  LOOP THRU THE STD LCS ONE BY ONE, DOING FITS ALONG THE WAY
       do 20 ilc=1,nstd

C   SUBROUTINE TO DO THE FITTING FOR THIS, THE ITH STDLC
C          call fitlc(nobs,p,v,ilc,nslc,sp,sv,xx,xsq)
          call fitlc(xx,xsq)
C   THIS PREVENTS NAN OUTPUT
          if(xsq.ne.xsq) xsq=78.9
C   WRITE COEFFICIENTS TO cOUTFILE
          write(10,333) ilc,xx(1),xx(2),xx(3),xsq,nslc
 333      format(i3,3f8.3,e12.4,i5)

C   SUBROUTINE TO MAKE FULL FITTED LIGHT CURVE FOR ITH STDLC
          call mklc(ilc,nslc,sp,sv,xx,spf,svf)
C   WRITE FITTED LIGHT CURVE TO aOUTFILE (END IT WITH (0,0))
C    THREE TIMES THRU: (1) FOR PHASE < 0, (2) FOR 0 < PHASE < 1,
C    (3) FOR PHASE > 1.
          do 37 klm=1,3
             do 38 j=1,nslc
       if(klm.eq.1.and.spf(j).ge.0.5) write(11,383) spf(j)-1.,svf(j)
       if(klm.eq.2.) write(11,383) spf(j),svf(j)
       if(klm.eq.3.and.spf(j).le.0.5) write(11,383) spf(j)+1.,svf(j)
 383            format(2f8.3)
 38             continue
 37          continue
             write(11,388)
 388         format(' 0. 0.')

 20       continue
C  PUT AWAY THE COEEFICIENTS AND LC-ALL FILES... DONE MAKING THEM         
 21    close(11)
       close(10)


C  FIND WHICH ONES HAD THE SMALLEST CHI-SQUARED... SORT & PICK 4 BEST
C   ILC = ITH LC IN ORDER INPUT
C   ITH(I) = INDEX AFTER SORTING (1-6 = BEST TO WORST)
C          write(6,100) outfile
       open(10,file='c'//outfile,status='old')
       do 40 i=1,nstd
          read(10,*) ilc,xx1(i),xx2(i),xx3(i),xxsq(i),nslc
          ith(ilc) = ilc
C          write(6,413) ilc,xxsq(ilc)
C 413      format('  Istd, chisq = ',i3,e12.4)
 40       continue
       close(10)

       call sort(nstd,xxsq,xx1,xx2,xx3,ith)

       write(6,415)
 415   format(' ')
       do 44 i=1,4
       write(6,414) ith(i),xx1(i),xx2(i),xx3(i),xxsq(i)
 414   format('Islc, A, P0, Vi, chisq = ',i3,3f8.3,e12.4)
 44       continue


C  READ THE FOUR BEST LCS OUT OF aOUTFILE AND WRITE THEM 
C   (IN ORDER) TO bOUTFILE  (PHI, V1, V2, V3, V4)
C          write(6,100) outfile
       open(11,file='a'//outfile,status='old')
       open(12,file='b'//outfile,status='new')
C    GET VMAX/VMIN WHILE WE ARE HERE
       vmax = -99.
       vmin = 99.
       do 60 j=1,4
          k=0
          n00 = 0
          do 50 i=1,32000
             read(11,*,end=51) x,y
             if(x.eq.0..and.y.eq.0.) then
               n00 = n00+1
C               write(6,444) i,n00,j,ith(j)
C 444           format(' i, n00, j, ith(j) = ',4i6)
               goto 50
             endif
             if(ith(j).eq.n00+1) then 
               k = k+1
               slcx(j,k) = x
               slcy(j,k) = y
             else
               if(ith(j).lt.n00+1) then
                 goto 51
               else
                 goto 50
               endif
             endif
C           if(k.le.4.or.k.ge.97) write(6,*) j,ith(j),i,k,slcx(j,k),slcy(j,k)
             vmax = max(vmax,slcy(j,k))
             vmin = min(vmin,slcy(j,k))             
 50          continue

 51       rewind(11)
 60       continue
       do 70 i=1,k
          write(12,200) slcx(1,i),slcy(1,i),slcy(2,i),slcy(3,i),slcy(4,i)
 200      format(5f8.3)
 70       continue
       close(12)
       close(11)


C  WRITE THE DATA INTO dOUTFILE (VMAG, PTYP, PHI1, PHI2, PHI3, PHI4, VERR)
C   IE ONE FOR EACH OF 4 BEST FITTED LCS -- PHI(FIT) = PHI(O) - DPHI
       open(17,file='d'//outfile,status='new')
       do 83 j=1,nobs
          pfit1 = p(j) + xx2(1)
          if(pfit1.lt.0.) pfit1 = pfit1+1.0
          if(pfit1.ge.1.) pfit1 = pfit1-1.0
          pfit2 = p(j) + xx2(2)
          if(pfit2.lt.0.) pfit2 = pfit2+1.0
          if(pfit2.ge.1.) pfit2 = pfit2-1.0
          pfit3 = p(j) + xx2(3)
          if(pfit3.lt.0.) pfit3 = pfit3+1.0
          if(pfit3.ge.1.) pfit3 = pfit3-1.0
          pfit4 = p(j) + xx2(4)
          if(pfit4.lt.0.) pfit4 = pfit4+1.0
          if(pfit4.ge.1.) pfit4 = pfit4-1.0
          write(17,222) v(j),nptp(j),pfit1,pfit2,pfit3,pfit4,ve(j)
 222      format(f7.3,i4,5f8.3)
 83       continue
       close(17)


C  MAKE A SUPERMONGO PLOT FILE AND THE OBSERVED DATA FILE FOR IT --
C    THIS ONE IS FOR THE 4V-FIT DIAGNOSTIC... USE TO PICK A BEST FIT.
 59    call mksm(nobs,nslc,vmax,vmin,xxsq,outfile)
       write(6,345)
 345   format(/,'SM: ')
       write(6,346) outfile
 346   format('  macro read s',a12)
       write(6,347)
 347   format('  lc')

C  CHOSE THE BEST LC BY EYE (DEF IS 1)
       write(6,678)
 678   format(/'Enter plot number (1-4) for final LC:  ',$)
       read(5,*) nfin
       if(nfin.lt.1.or.nfin.gt.4) nfin = 1

C  WRITE FINAL OBSERVED (PHASE-SHIFTED) DATA FILE (PHI V NPTP),
C   AND COMPUTE LC PARAMETERS BASED ON OBSERVED DATA.
       open(33,file=origf,status='old')
       open(17,file='g'//outfile,status='new')
C   WHILE WERE HERE, COMPUTE THE BRIGHTEST I POINT AND FAINTEST V
C    OF THE OBSERVED DATA FOR FINAL SM PLOT (SEE LAST SECTION OF MAIN PROG).
C    AND ALSO WORK ON COMPUTATION OF COLOR AT (PHASE-SHIFTED) MIN LIGHT.
       nallvi = 0
       ybrt =  88.
       yfnt = -99.
       hjdmax = -9999.
       hjd0 = -9999.
       cmin = 0.
       csq = 0.
       nmin = 0
C       read(33,*) origper
C       do 183 j=1,nobs
       do 183 j=1,99999
C	  read(33,*,end=184) pold,vold,nold,hjdold,veold,bold,beold
C     THIS FORMAT IS THE SAME AS USED BY stdlc_per10_export.f
          read(33,*,end=184) hjdold, nold, vold, veold, bold, beold
          nallvi = nallvi + 1
C     IN THIS FORMAT, WE NEED TO COMPUTE A PROVISIONAL PHASE, PI,
C       BASED ON THE FIRST-ENCOUNTERED OBSERVATION (==> PI=0).
          if(j.eq.1) hjd0 = hjdold
          cyc = (hjdold-hjd0)/period
          pold = cyc - float( int(cyc) )
          if(pold.lt.0.0) pold = pold+1.0
          if(pold.ge.1.0) pold = pold-1.0
C     BULLETPROOF AGAINST OOB DATA IN INPUT FILE
          if(vold.gt.yfnt .and. vold.lt.40.0) yfnt = vold
          if(bold.gt.yfnt .and. bold.lt.40.0) yfnt = bold
C          yfnt = max(yfnt,bold,vold)
C          if(bold.gt.40.0) goto 1033
C          if(vold.gt.40.0) goto 1033
          ybrt = min(ybrt,bold,vold)
C          if(bold.lt.40.0) then
C            if((bold+offi).lt.ybrt) ybrt = bold+offi
C           else
C            ybrt = min(ybrt,vold)
C          endif
C          if(vold.lt.40.0) then
C            if((vold+offi).lt.ybrt) ybrt = vold+offi
C           else
C            ybrt = min(ybrt,bold)
C          endif
 1033     pold = pold + xx2(nfin)
          if(pold.lt.0.0) pold = pold+1.0
          if(pold.ge.1.0) pold = pold-1.0
          if(hjdold.gt.hjdmax) lasthjd = j
C   COLOR AT MIN LIGHT
          if(pold.lt.0.5 .or. pold.gt.0.8) goto 645
          if(vold.gt.filtmag1 .or. vold.lt.filtmag2) goto 645
          if(bold.gt.filtmag1 .or. bold.lt.filtmag2) goto 645
          cdat = vold - bold
          cmin = cmin + cdat
          csq = csq + cdat*cdat
          nmin = nmin +1
C   APPLY ERRORBAR FILTERS
 645      if(vold.gt.filtmag1 .or. vold.lt.filtmag2) veold = filterr
          if(bold.gt.filtmag1 .or. bold.lt.filtmag2) beold = filterr
          write(17,122) hjdold,pold,vold,veold,bold,beold,nold,outfile
 122      format(f14.5,5f7.3,i4,1x,a8)
C          p(j) = p(j) + xx2(nfin)
C          if(p(j).lt.0.) p(j) = p(j)+1.0
C          if(p(j).ge.1.) p(j) = p(j)-1.0
C          if(hjd(j).gt.hjdmax) lasthjd = j
CC          write(17,122) p(j),v(j),nptp(j),hjd(j),verr(j)
C          write(17,122) hjd(j),p(j),v(j),verr(j),nptp(j),outfile
C 122      format(f14.5,3f7.3,i4,1x,a8)
 183      continue
 184    close(17)

C  HERE IS WHERE WE FINALIZE THE LC PARAMS FROM THE OBSERVED DATA POINTS
C   COMPUTE COLOR AT MINIMUM LIGHT AND STATS.  OMIT BAD DATA PTS (EG NO I OBS)
        if(nmin.lt.1) then
          cmino = 99.999
          emino = 88.888
          goto 1435
        endif
        cmino = cmin/float(nmin)
        if(nmin.eq.1) then
          emino = 99.999
          goto 1435
        endif
        sd = sqrt( (csq - cmin*cmin/float(nmin)) / (float(nmin-1)) )
        emino = sd/sqrt(float(nmin))
C  WRITE THE DERIVED DATA TO THE DATA OUTPUT FILE
 1435  open(17,file='rd.'//datfile)
       do 147 i=1,32000
          read(17,100,end=148) infile
 147       continue
 148   write(6,4555)
 4555  format(/,'Light Curve Parameters from Observed Data Points:')
       write(6,455)
 455   format('outfile    <V>o   Vmaxo  Vmino   Amplo  dPro      Emaxo     Nobs  <I>o  V-Imin  sem   Nmin')
       write(6,456)  outfile,vimo,vmaxo,vmino,amplo,dphiro,emaxo,nobs,iimo,cmino,emino,nmin
       write(17,456) outfile,vimo,vmaxo,vmino,amplo,dphiro,emaxo,nobs,iimo,cmino,emino,nmin
 456   format(a10,5f7.3,f14.5,i5,3f7.3,i4)
       close(17)


C  WRITE FINAL FITTED DATA FILE, fOUTFILE (phi V),
C   AND COMPUTE LC PARAMETERS BASED ON FITTED DATA.
       open(17,file='f'//outfile,status='new')
       vminf = -99.
       vmaxf = 99.
       k = 0
       do 187 j=1,2*nslc+1
          write(17,127) slcx(nfin,j),slcy(nfin,j)
 127      format(2f8.3)
          if(slcx(nfin,j).lt.0..or.slcx(nfin,j).ge.1.) goto 187
          k = k+1
          spf(k) = slcx(nfin,j)
          svf(k) = slcy(nfin,j)
C   NOTE: VMIN/MAX IN MAG/LC SENSE (SMALL NUMBERS IS MAX).
          if(slcy(nfin,j).gt.vminf) then
            vminf = slcy(nfin,j)
            jmin = j
          endif
          if(slcy(nfin,j).lt.vmaxf) then
            vmaxf = slcy(nfin,j)
            jmax = j
          endif
 187      continue
       amplf = vminf-vmaxf
       dphirf = slcx(nfin,jmax)-slcx(nfin,jmin)
       if(dphirf.lt.0.) dphirf = dphirf+1.
       if(dphirf.ge.1.) dphirf = dphirf-1.
       call getimean(nslc,spf,svf,vimf)
       emaxf = hjd(lasthjd) - p(lasthjd)*period
       close(17)

C  WRITE THE DERIVED LC PARAMS FROM THE FIT TO THE OUTPUT FILE
       open(17,file='fd.'//datfile)
       do 247 i=1,32000
          read(17,100,end=248) infile
 247      continue
 248   rms = sqrt(xxsq(nfin))
       write(6,4556)
 4556  format(/,'Light Curve Parameters from FIT:')
       write(6,2455) 
 2455  format('infile     <V>f   Vmaxf  Vminf   Amplf  dPrf      Emaxf   Nobs slc RMS    Pshft')
       write(6,2456) outfile,vimf,vmaxf,vminf,amplf,dphirf,emaxf,nobs,ith(nfin),rms,xx2(nfin)
       write(17,2456) outfile,vimf,vmaxf,vminf,amplf,dphirf,emaxf,nobs,ith(nfin),rms,xx2(nfin)
 2456   format(a10,5f7.3,f14.5,2i3,2f7.3)
       close(17)

C  MAKE A SUPERMONGO PLOT FILE AND THE OBSERVED DATA FILE FOR IT --
C    THIS ONE IS FOR THE **FINAL** V AND I LIGHT CURVES!!
       open(17,file='z'//outfile)
       call mksmvi(nallvi,nslc,nfin,ybrt,yfnt,offi,outfile)
       write(6,3345)
 3345   format(/,'SM: ')
       write(6,3346) outfile
 3346   format('  macro read z',a12)
       write(6,3347)
 3347   format('  lc')
       close(17)


       goto 1
 2     stop 
       end


C------------------------------------------------------------------------
C  THIS BIT TAKES THE STD LC AND THE FITTED COEFFICIENTS AND
C   CONSTRUCTS A FITTED LIGHT CURVE
      subroutine mklc(ilc,nslc,sp,sv,xx,spf,svf)

      real sp(256),sv(10,256),spf(256),svf(256)
      real xx(3)

      do 10 j=1,nslc
C   PHASE IS PHASE...NB WE MAKE THE FITTED LC IN PLACE, AND WILL SHIFT
C    THE OBSD LC TO MATCH THE PHASE (SINCE THE STD LC HAS A
C    BETTER-DEFINED MAXIMUM)
         spf(j) = sp(j)
         svf(j) = xx(1)*sv(ilc,j)+xx(3)
 10      continue
      return
      end


C------------------------------------------------------------------------
C  THIS SECTION FITS THE STD LC TO THE OBSERVED LC
C   
C      subroutine fitlc(nobs,p,v,ilc,nslc,sp,sv,xx,xsq)
      subroutine fitlc(xx,xsq)

      real p(256),v(256),sp(256),sv(10,256),ve(256)
      real pp(4,3),yy(4),xx(3),imag(256),ierr(256)
      character*1 ckey
      common/lcdat1/nobs,p,v,ve,imag,ierr
      common/lcdat2/nslc,ilc,sp,sv
      
C  IVERB=1 GIVES REAL TIME SCREEN OUTPUT OF FITTING, IVERB=0 DOESNT.
      iverb=0

C  SET NDIM, MP=4 AND NP=3, ESSENTIALLY SOFTWIRING THE SIMPLEX CODE
C   'AMOEBA' TO THESE VALUES.
      ndim = 3
      mp = ndim+1
      np = ndim
      iter = 0

C  START BY GETTING INITIAL GUESSES FOR SIMPLEX
      vmax = -99.
      vmin = 99.
      do 10 i=1,nobs
         if(v(i).lt.vmin) then
           vmin = v(i)
           pmin = p(i)
         endif
         if(v(i).gt.vmax) vmax = v(i)
 10      continue
C   BEST GUESS IS A=VMAX-VMIN, P0=P(VMIN), VI=(VMAX+VMIN)/2.
C   SET PP(I,J) TO BEST GUESS
      do 20 i=1,4
         pp(i,1) = vmax-vmin
         pp(i,2) = pmin
         pp(i,3) = (vmax+vmin)/2.
 20      continue
C   AND DO OTHER 3 AMOEBA GUESSES AS 1-D DEVIATES FROM THIS.
      pp(2,1) = pp(2,1)+0.2
      pp(3,2) = pp(3,2)+0.2
      pp(4,3) = pp(4,3)+0.4

C  REALITY CHECK
      if(iverb.eq.1) then
      do 22 i=1,4
         write(6,220) i,pp(i,1),pp(i,2),pp(i,3)
 220     format('PP(',i1,',N) = ',3e12.4)
 22      continue
      endif

C  NOW EVALUATE THE FUNCTION TO BE MINIMIZED FOR EACH OF THE 4
C   GUESSES IN PP, AND CALL THE VECTOR YY.
      do 32 i=1,4
         do 31 j=1,3
            xx(j) = pp(i,j)
 31         continue
         yy(i) = funk(xx)
 32      continue

C  REALITY CHECK
      if(iverb.eq.1) then
      do 37 i=1,4
         write(6,370) i,yy(i)
 370     format('YY(',i1,') = ',e12.4)
 37      continue

C  HOLD UP HERE SO I CAN LOOK AT PP AND YY IN REAL TIME
      write(6,388)
 388  format('Hit any key to continue:  ',$)
      read(5,111) ckey
 111  format(a)
      endif

C  SET THE FRACTIONAL CONVERGENCE TOLERANCE (IN UNITS OF FUNCTION VALUE)
C      ftol = 0.00001
C   SOFT VERSION!
      ftol = 0.0001

C  CALL THE SIMPLEX OPTIMIZATION ROUTINE AMOEBA
      call amoeba(pp,yy,mp,np,ndim,ftol,funk,iter)

C  RESULTS OF OPTIMIZATION ARE 4 VECTORS REAL CLOSE TO THE SAME
      suma = 0.
      sump0 = 0.
      sumvi = 0.
      if(iverb.eq.1) write(6,320) 
 320  format('          N       a        p0       Vi')
      do 40 i=1,4
         if(iverb.eq.1) write(6,330) i,pp(i,1),pp(i,2),pp(i,3)
 330     format(10x,i1,3f8.3)
         suma = suma+pp(i,1)
         sump0 = sump0+pp(i,2)
         sumvi = sumvi+pp(i,3)
 40      continue
C  FOR NOW, TAKE THE MEAN OF THE 4 VECTORS AS THE SOLUTION.
      xx(1) = suma/4.
      xx(2) = sump0/4.
      xx(3) = sumvi/4.
C  EVALUATE IT AT THE FINAL SOLUTION, FOR LATER CHI-SQUARED EVALUATION.
      xsq = funk(xx)

C  FINAL VECTOR (THIS HAS *NOT* BEEN RESTARTED NEAR FINAL SOLUTION)
      if(iverb.eq.1) write(6,420) xx(1),xx(2),xx(3),xsq,iter
 420  format(/,'First Soln: ',3f8.3,' chisq: ',e12.4,
     & ' Niter: ',i5)

C  NOW RESTART THE SOLUTION USING AN INITIAL GUESS HAVING FIRST
C   ROW SAME (IE, FIRST ROW OF BEST SOLUTION) AND ROWS 2-4 HAVING
C   1.1*BEST SOLUTION, WITH A LITTLE SHUFFLING OF COEFFS.

C  HERES THE SHUFFLE
      pp(2,1) = pp(3,1)
      pp(3,2) = pp(4,2)
      pp(4,3) = pp(2,3)
C  AND HERE'S THE 1.1*
      do 55 i=2,4
         do 54 j=1,3
            pp(i,j) = 1.1*pp(i,j)
 54         continue
 55      continue

C  NOW EVALUATE THE FUNCTION TO BE MINIMIZED FOR EACH OF THE 4
C   GUESSES IN PP, AND CALL THE VECTOR YY.
      do 62 i=1,4
         do 61 j=1,3 
            xx(j) = pp(i,j)
 61         continue
         yy(i) = funk(xx)
 62      continue

C  CALL THE SIMPLEX OPTIMIZATION ROUTINE AMOEBA
      iter = 0
      call amoeba(pp,yy,mp,np,ndim,ftol,funk,iter)

C  RESULTS OF OPTIMIZATION ARE 4 VECTORS REAL CLOSE TO THE SAME
      suma = 0.
      sump0 = 0.
      sumvi = 0.
      if(iverb.eq.1) write(6,320) 
      do 70 i=1,4
         if(iverb.eq.1) write(6,330) i,pp(i,1),pp(i,2),pp(i,3)
         suma = suma+pp(i,1)
         sump0 = sump0+pp(i,2)
         sumvi = sumvi+pp(i,3)
 70      continue
C  FOR NOW, TAKE THE MEAN OF THE 4 VECTORS AS THE SOLUTION.
      xx(1) = suma/4.
      xx(2) = sump0/4.
      xx(3) = sumvi/4.
C  EVALUATE IT AT THE FINAL SOLUTION, FOR LATER CHI-SQUARED EVALUATION.
      xsq = funk(xx)
C  IF MAX NUMBER OF ITERATIONS IS EXCEEDED IN AMOEBA (IE IT DIDN'T CONVERGE),
C   USE THESE (BEST GUESS) VALUES OF XX(1-3) BUT SET CHISQ=99.99 SO THIS
C   SOLUTION DOESN'T GET PICKED BY THE AUTOMATIC ROUTINES.
	if(iter.gt.4999) xsq = 99.99

C  FINAL VECTOR (THIS *HAS* BEEN RESTARTED NEAR FINAL SOLUTION)
      write(6,720) xx(1),xx(2),xx(3),xsq,iter,ilc
 720  format('Final Soln: ',3f8.3,' chisq: ',e12.4,
     & ' Niter: ',i5,' Ilc: ',i2)

      return
      end


C----------------------------------------------------------------------
C  FUNCTION TO EVALUATE THE CHI-SQUARED RESIDUALS BETWEEN THE
C   FITTED AND OBSERVED LCS GIVEN THE 3-ELEMENT VECTOR (A,P0,VI)
      function funk(xx)

      real p(256),v(256),sp(256),sv(10,256),ve(256)
      real xx(3),imag(256),ierr(256)
C      common/lcdat/nobs,p,v,nslc,sp,sv,ilc,ve
      common/lcdat1/nobs,p,v,ve,imag,ierr
      common/lcdat2/nslc,ilc,sp,sv

C   NO NEED TO WRAP, SINCE WE HAVE A SLC PT AT (P,V) = (1.0,V(1))
      sresq = 0.
      swtsq = 0.
C   FOR EACH POINT IN THE OBSERVED SPECTRUM, COMPUTE THE RESIDUAL
C    = DIFFERENCE FROM EXPECTED VALUE AT THAT PHASE GIVEN PRESENT
C    BEST GUESS AT COEFFICIENTS A,P0,VI
C       WRITE(6,777) NOBS,NSLC,OP
C 777   FORMAT('Nobs, Nslc = ',2I5)
      do 10 i=1,nobs
C    REMEMBER, THE OBSD LC PHASE PTS ARE SHIFTED BY P0=XX(2) HERE.
         op = p(i)+xx(2)
         if(op.lt.0.) op=op+1.
         if(op.gt.1.) op=op-1.
         do 9 j=2,nslc
C   STEP THROUGH STD LC POINTS TILL PHASE OF OBSD PT IS BTWN 2 STD PTS
            if(op.lt.sp(j).and.op.ge.sp(j-1)) then
              svh = xx(1)*sv(ilc,j)+xx(3)
              svl = xx(1)*sv(ilc,j-1)+xx(3)
C   LINEAR INTERPOLATION AT OP TO GET VSTD(OP) AT CURRENT BEST FIT CURVE
              vint = vintrp(op,sp(j),sp(j-1),svh,svl)
              res = v(i) - vint
C   COMUPUTE SUM OF SQUARED RESIDUALS... THAT'S WHAT WE ARE MINIMIZING
C	SEVERAL WEIGHTING OPTIONS: NONE, 1/ERR, 1/ERR^2
C              sresq = sresq + res*res
C              sresq = sresq + res*res/ve(i)
C              swtsq = swtsq + 1.0/ve(i)
              sresq = sresq + res*res/(ve(i)*ve(i))
              swtsq = swtsq + 1.0/(ve(i)**2)

C       WRITE(6,778) i,p(i),v(i),j,op,vint,res
C 778   FORMAT('i,P(i),V(i),j,OP,Vint,resid = ',i3,2f8.3,i5,3f8.3)
            endif
 9          continue
 10      continue
C   MAKE IT LIKE CHI-SQUARED BY DIVIDING BY NUMBER OF OBSERVED POINTS
C      funk = sresq/float(nobs)
C   MAKE IT LIKE CHI-SQUARED BY DIVIDING BY SUM OF WTS
      funk = sresq/swtsq
C      write(6,222) funk
C 222  format('FUNK = ',e12.4)
      return
      end


C-----------------------------------------------------------------------
C  FUNCTION TO INTERPOLATE BETWEEN 2 STD LC PTS
      function vintrp(p,sph,spl,svh,svl)

      dx = sph-p
      xxx = sph-spl
      yyy = svh-svl
      vintrp = svh-dx*yyy/xxx
      return
      end

C------------------------------------------------------------------------
C  THIS IS THE SIMPLEX ROUTINE AMOEBA FROM 
C   /home/physun/recipes/recipes_f/recipes/amoeba.f = P292 OF NUM REC.
      SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      INTEGER iter,mp,ndim,np,NMAX,ITMAX
      REAL ftol,p(mp,np),y(mp),funk
C      PARAMETER (NMAX=20,ITMAX=50000)
      PARAMETER (NMAX=20,ITMAX=5000)
      EXTERNAL funk
CU    USES amotry,funk
      INTEGER i,ihi,ilo,inhi,j,m,n
      REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry
      iter=0
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue
        return
      endif
C      if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba'
      if(iter.ge.ITMAX) then
        write(6,789)
 789    format('WARNING: ITMAX exceeded in amoeba, USING PRESENT VALUES')
        return
      endif
      iter=iter+2
      ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0)
      if (ytry.le.y(ilo)) then
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              y(i)=funk(psum)
            endif
16        continue
          iter=iter+ndim
          goto 1
        endif
      else
        iter=iter-1
      endif
      goto 2
      END


C------------------------------------------------------------------------
C  A FUNCTION USED BY AMOEBA
      FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac)
      INTEGER ihi,mp,ndim,np,NMAX
      REAL amotry,fac,p(mp,np),psum(np),y(mp),funk
      PARAMETER (NMAX=20)
      EXTERNAL funk
CU    USES funk
      INTEGER j
      REAL fac1,fac2,ytry,ptry(NMAX)
      fac1=(1.-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      ytry=funk(ptry)
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)
12      continue
      endif
      amotry=ytry
      return
      END


C------------------------------------------------------------------------
C  MAKES A SUPERMONGO FILE MACRO FOR POTTING ALL THE LCS ON ON GRAPH
      subroutine mksm(nobs,nslc,vmax,vmin,xxsq,outfile)

      character*50 outfile
      dimension xxsq(4)
     
      open(17,file='s'//outfile,status='new')

C  READ DATA INTO SM
C   FIRST THE REAL STAR DATA (V NPT P1 P2 P3 P4 VERR)
      write(17,190)
 190  format('lc     erase')
      write(17,200) outfile
 200  format('       data d',a50)
      write(17,210) nobs
 210  format('       lines 1 ',i6)
      write(17,220) 
 220  format('       read { yo 1  ptp 2  xo1 3  xo2 4  xo3 5  xo4 6  eyo 7 }')
C   NOW THE 4 BEST STDLCS (PHI V1 V2 V3 V4)
      write(17,222) outfile
 222  format('       data b',a50)
      write(17,224) 2*nslc+1
 224  format('       lines 1 ',i6)
      write(17,226) 
 226  format('       read { xf 1  yf1 2  yf2 3  yf3 4  yf4 5 }')

C  SET UP WINDOW FOR FIRST (UL) FIT
      write(17,230) vmax+0.2,vmin-0.1
 230  format('       limits -0.49 1.49 ',2f8.3)
C      write(17,240) 
C 240  format('       window 1 2 2 2')
      write(17,240) 
 240  format('       ticksize  0.05 0.5   0.05 0.5')
      write(17,250) 
 250  format('       location 3000 17000 16000 29000')
      write(17,260) 
 260  format('       expand 1.4')
      write(17,270) 
 270  format('       box 0 1')
      write(17,280) 
 280  format('       ptype ptp')
      write(17,283) 
 283  format('       points xo1 yo')
      write(17,284) 
 284  format('       errorbar xo1 yo eyo 2')
      write(17,285) 
 285  format('       errorbar xo1 yo eyo 4')
      write(17,286)
 286  format('       connect xf yf1')
C      write(17,293) 
C 293  format('       xlabel \phi')
      write(17,296) 
 296  format('       ylabel V')
      write(17,297) 
 297  format('       relocate (3500 16500)')
      write(17,298) 
 298  format('       expand 1.1')
      write(17,299) sqrt(xxsq(1))
 299  format('       putlabel 9 "1: rms = ',f8.4,'"')

C  SET UP WINDOW FOR SECOND (UR) FIT
      write(17,330) vmax+0.2,vmin-0.1
 330  format('       limits -0.49 1.49 ',2f8.3)
C      write(17,340) 
C 340  format('       window 2 2 2 2')
      write(17,240) 
      write(17,350) 
 350  format('       location 17000 31000 16000 29000')
      write(17,260) 
      write(17,370) 
 370  format('       box 0 0')
      write(17,380) 
 380  format('       ptype ptp')
      write(17,383) 
 383  format('       points xo2 yo')
      write(17,384) 
 384  format('       errorbar xo2 yo eyo 2')
      write(17,385) 
 385  format('       errorbar xo2 yo eyo 4')
      write(17,386)
 386  format('       connect xf yf2')
C      write(17,393) 
C 393  format('       xlabel \phi')
C      write(17,396) 
C 396  format('       ylabel V')
      write(17,397) 
 397  format('       relocate (17500 16500)')
      write(17,298) 
      write(17,399) sqrt(xxsq(2))
 399  format('       putlabel 9 "2: rms = ',f8.4,'"')

C  SET UP WINDOW FOR THIRD (LL) FIT
      write(17,430) vmax+0.2,vmin-0.1
 430  format('       limits -0.49 1.49 ',2f8.3)
C      write(17,440) 
C 440  format('       window 1 1 2 2')
      write(17,240) 
      write(17,450) 
 450  format('       location 3000 17000 3000 16000')
      write(17,260) 
      write(17,470) 
 470  format('       box 1 1')
      write(17,480) 
 480  format('       ptype ptp')
      write(17,483) 
 483  format('       points xo3 yo')
      write(17,484) 
 484  format('       errorbar xo3 yo eyo 2')
      write(17,485) 
 485  format('       errorbar xo3 yo eyo 4')
      write(17,486)
 486  format('       connect xf yf3')
      write(17,493) 
 493  format('       xlabel \phi')
      write(17,496) 
 496  format('       ylabel V')
      write(17,497) 
 497  format('       relocate (3500 3500)')
      write(17,298)
      write(17,499) sqrt(xxsq(3))
 499  format('       putlabel 9 "3: rms = ',f8.4,'"')

C  SET UP WINDOW FOR FOURTH (LR) FIT
      write(17,530) vmax+0.2,vmin-0.1
 530  format('       limits -0.49 1.49 ',2f8.3)
C      write(17,540) 
C 540  format('       window 2 1 2 2')
      write(17,240) 
      write(17,550) 
 550  format('       location 17000 31000 3000 16000')
      write(17,260) 
      write(17,570) 
 570  format('       box 1 0')
      write(17,580) 
 580  format('       ptype ptp')
      write(17,583) 
 583  format('       points xo4 yo')
      write(17,584) 
 584  format('       errorbar xo4 yo eyo 2')
      write(17,585) 
 585  format('       errorbar xo4 yo eyo 4')
      write(17,586)
 586  format('       connect xf yf4')
      write(17,593) 
 593  format('       xlabel \phi')
C      write(17,596) 
C 596  format('       ylabel V')
      write(17,597) 
 597  format('       relocate (17500 3500) ')
      write(17,298)
      write(17,599) sqrt(xxsq(4))
 599  format('       putlabel 9 "4: rms = ',f8.4,'"')

      write(17,260)
      write(17,600) 
 600  format('       relocate (17000 30000)')
      write(17,605) outfile
 605  format('       putlabel 5 ',a12)

      close(17)
      return
      end




C------------------------------------------------------------------------
C  MAKES A SUPERMONGO FILE MACRO FOR POTTING ALL THE LCS ON ON GRAPH
      subroutine mksmvi(nobs,nslc,nfin,ybrt,yfnt,offi,outfile)

      character*50 outfile
C      dimension xxsq(4)
     
C      open(17,file='s'//outfile,status='new')

C  READ DATA INTO SM
C   FIRST THE REAL STAR DATA (V NPT P1 P2 P3 P4 VERR), SHIFTED FOR FITPHASE
      write(17,190)
 190  format('lc     erase')
      write(17,200) outfile
 200  format('       data g',a50)
      write(17,210) nobs
 210  format('       lines 1 ',i6)
      write(17,220) 
 220  format('       read { phi 2  vmag 3  verr 4  imag 5  ierr 6  ptp 7 }')
C 220  format('       read { yo 1  ptp 2  xo1 3  xo2 4  xo3 5  xo4 6  eyo 7 }')
C   NOW THE FINAL FITTED STDLCS (PHI V)
      write(17,222) outfile
 222  format('       data f',a50)
      write(17,224) 2*nslc+1
 224  format('       lines 1 ',i6)
      write(17,226) 
 226  format('       read { phif 1  vf 2 }')
C  SET UP WINDOW
      write(17,230) yfnt+0.2,ybrt-0.2
 230  format('       limits -0.3 1.3 ',2f8.3)
      write(17,250) 
 250  format('       location 5000 30000 4000 30000')
      write(17,253) 
 253  format('       ticksize  0.1 0.5  0.1 0.5')
      write(17,260) 
 260  format('       expand 1.8')
      write(17,270) 
 270  format('       box 1 1')
      write(17,280) 
 280  format('       ptype ptp')
C  WRITE OBSD VMAG POINTS IN BLUE
      write(17,283) 
 283  format('       ctype blue')
      write(17,2283) 
2283  format('       points phi vmag')
      write(17,284) 
 284  format('       errorbar phi vmag verr 2')
      write(17,285) 
 285  format('       errorbar phi vmag verr 4')
      write(17,290)
 290  format('       relocate (27500 10000)')
      write(17,293)
 293  format('       putlabel 5 V')
C  WRITE OBSD IMAG POINTS IN RED
C      write(17,380) offi 
C 380  format('       set imag = imag + ',f6.3)
      write(17,3383) 
3383  format('       ctype red')
      write(17,383) 
 383  format('       points phi imag')
      write(17,384) 
 384  format('       errorbar phi imag ierr 2')
      write(17,385) 
 385  format('       errorbar phi imag ierr 4')
      write(17,390)
 390  format('       relocate (27500 25000)')
        write(17,393) 
 393    format('       putlabel 5 I')
C      if(offi.ge.0.0) then
C        write(17,393) offi
C 393    format('       putlabel 5 I+',f4.2)
C      else
C        write(17,394) offi
C 394    format('       putlabel 5 I',f5.2)
C      endif
C  WRITE FITTED VMAG CURVE IN GREEN
      write(17,483) 
 483  format('       ctype green')
      write(17,4483) 
4483  format('       connect phif vf')
      write(17,490)
 490  format('       relocate (6600 11000)')
      write(17,493)
 493  format('       putlabel 5 fit')
C  FINISH UP
      write(17,586)
 586  format('       ctype black')
      write(17,593) 
 593  format('       xlabel Phase')
      write(17,596) 
 596  format('       angle 90')
      write(17,597) 
 597  format('       relocate (1500 17000)')
      write(17,598) 
 598  format('       putlabel 5 mag')
      write(17,5598) 
5598  format('       angle 0')
      write(17,599) 
 599  format('       relocate (17500 31000)')
      write(17,5599) outfile
5599  format('       putlabel 9 ',a8)

C      close(17)
      return
      end



C  ROUTINE FROM "NUMERICAL RECIPES" P.231.  SORTS A KEY ARRAY (RA) AND 
C  4 COMPANION ARRAYS (CA & ICA) INTO ASCENDING ORDER OF THE KEY USING 
C  THE "HEAPSORT" AGORITHM.			AL 3/89
       subroutine sort(n,ra,ca,cb,cc,ica)

       dimension ra(n),ca(n),cb(n),cc(n)
       integer ica(n)
       l = n/2+1
       ir = n
10     continue
          if(l.gt.1) then
            l = l-1
            rra = ra(l)
            rca = ca(l)
            rcb = cb(l)
            rcc = cc(l)
            irca = ica(l)
          else
            rra = ra(ir)
            rca = ca(ir)
            rcb = cb(ir)
            rcc = cc(ir)
            irca = ica(ir)
            ra(ir) = ra(1)
            ca(ir) = ca(1)
            cb(ir) = cb(1)
            cc(ir) = cc(1)
            ica(ir) = ica(1)
            ir = ir-1
            if(ir.eq.1) then
              ra(1) = rra
              ca(1) = rca
              cb(1) = rcb
              cc(1) = rcc
              ica(1) = irca
              return
            endif
          endif
          i = l
          j = l+l
20        if(j.le.ir) then
            if(j.lt.ir) then
              if(ra(j).lt.ra(j+1)) j = j+1
            endif
            if(rra.lt.ra(j)) then
              ra(i) = ra(j)
              ca(i) = ca(j)
              cb(i) = cb(j)
              cc(i) = cc(j)
              ica(i) = ica(j)
              i = j
              j = j+j
            else
              j = ir+1
            endif
          goto 20
          endif
          ra(i) = rra
          ca(i) = rca
          cb(i) = rcb
          cc(i) = rcc
          ica(i) = irca
       goto 10
       end



       subroutine getimean(nobs,phi,v,vimean)
C  THIS BEAST GETS THE INTENSITY-MEAN MAG OF THE STAR.

       dimension phi(256),v(256),x(256),p(256)
       co = 100000.

C  SET UP VALUES FOR I=1,N OBS.  PROTECT AGAINST BAD DATA VALS.
       nu = 0
       do 10 i=1,nobs
          if(v(i).gt.35.0.or.v(i).lt.4.0) goto 10
          nu = nu+1
          x(nu) = co * 10.**(-0.4*v(i))
          p(nu) = phi(i)
 10       continue

C     BULLET PROOF AGAINST LOW NUMBER OF OBS
       if(nu.lt.1) then
         vimean = 99.999
         goto 999
       endif
       if(nu.le.1) then
         vimean = 88.888
         goto 999
       endif

       call sort2(nu,p,x)

C  SET UP VALUES FOR THE N+1TH OBS = 1ST OBS AT PHI=PHI+1
       p(nu+1) = p(1)+1.0
       x(nu+1) = x(1)

       asum = 0.
       psum = 0.
       do 20 i=1,nu
          xa = (x(i+1)+x(i))/2.
          dp = p(i+1)-p(i)
          asum = asum + xa*dp
          psum = psum + dp
C        WRITE(6,*) I,X(I),XA,P(I),DP
 20       continue
C        WRITE(6,*) ASUM,PSUM
       a = asum/psum
       vimean = -2.5 * log10(a/co)
 999   return
       end

C  ************************* SORT2.FOR **************************
C  ROUTINE FROM "NUMERICAL RECIPES" P.231.  SORTS A KEY ARRAY (RA) AND 
C  2 COMPANION ARRAYS (CA & CB) INTO ASCENDING ORDER OF THE KEY USING 
C  THE "HEAPSORT" AGORITHM.			AL 3/89
       subroutine sort2(n,ra,ca)

       real ra(n),rra,ca(n),rca

       l = n/2+1
       ir = n
10     continue
          if(l.gt.1) then
            l = l-1
            rra = ra(l)
            rca = ca(l)
C            rcb = cb(l)
C            rcc = cc(l)
          else
            rra = ra(ir)
            rca = ca(ir)
C            rcb = cb(ir)
C            rcc = cc(ir)
            ra(ir) = ra(1)
            ca(ir) = ca(1)
C            cb(ir) = cb(1)
C            cc(ir) = cc(1)
            ir = ir-1
            if(ir.eq.1) then
              ra(1) = rra
              ca(1) = rca
C              cb(1) = rcb
C              cc(1) = rcc
              return
            endif
          endif
          i = l
          j = l+l
20        if(j.le.ir) then
            if(j.lt.ir) then
              if(ra(j).lt.ra(j+1)) j = j+1
            endif
            if(rra.lt.ra(j)) then
              ra(i) = ra(j)
              ca(i) = ca(j)
C              cb(i) = cb(j)
C              cc(i) = cc(j)
              i = j
              j = j+j
            else
              j = ir+1
            endif
          goto 20
          endif
          ra(i) = rra
          ca(i) = rca
C          cb(i) = rcb
C          cc(i) = rcc
       goto 10
       end
