CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  STDLC_PER.F -- FIND PERIOD OF AN RRLYR BY FOLDING THE 
C     TIME-MAG DATA BY A SERIES OF PERIODS, AND FOR EACH 
C     PERIOD, FITTING RR LC TEMPLATES.  WRITE OUT CHISQ
C     FOR EACH FIT -- MINIMA ARE CANIDATE PERIODS TO FOLLOW
C     UP WITH PFOLD.  COPIES AND ADAPTED FROM
C     /n/Kracken2/layden/IHALO/LCS/MONTEC/stdlc_mcr2.f
C            98MAR13(Friday!) ACL
C
C     THIS VERSION COPY/ADAPTED FROM
C      /n/Kracken5/layden/N6441/DAO/ASSEM/VARS/PFOLD/stdlc_per.f
C      FOR USE WITH 5 TEMPLATES, VONLY, AND HJD-GROUP SELECTION
C
C     THIS VERSION ADAPTED FROM 
C     /n/Kracken5/bbowes/M54/RRLYR/PERIOD/stdlc_per.f 
C     TO DO *WEIGHTED*-CHI-SQ FITTING, MORE ROBUST AGAINST PTS W/
C     LARGE EBARS.           ACL 98JUN15
C
C  EXPORT VERSION (19 MAY 2000)
C     COPIED FROM /usr/Kracken4/layden/M54/RRL/BEN/PERIOD/stdlc_per10.f
C     AND ADAPTED FOR GENERIC RUN DATES INPUTED FROM FILE (stdlc_per10.runs).
C
C  IF YOU USE THIS PROGRAM IN ANY PUBLISHED MATERIAL, PLEASE EMAIL
C     EMAIL ME (layden@baade.bgsu.edu) AND IF POSSIBLE, SEND ME A COPY
C     OF THE RESULTING PUBLICATION.  I WOULD BE PLEASED TO HEAR OF ANY
C     COMMENTS YOU HAVE REGARDING THE PROGRAM (THOUGH KEEP IN MIND I
C     AM NOT A PROFESSIONAL PROGRAMMER, AND THIS PROGRAM WAS NEVER REALLY
C     INTENDED FOR PUBLIC USE!).
C                                        ANDY LAYDEN, 19 MAY 2000
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       character infile*50, slcfile*50, outfile*50
C  GOOD FOR UP TO 256 POINTS IN OBSD AND STD LCS, UP TO 10 STDLCS
       real hjd(256),vmag(256),imag(256),v(256),p(256),verr(256),ve(256)
       real sv(10,256),sp(256),xx(3)
       parameter (maxruns=32)
       real tbeg(maxruns), tend(maxruns)
       integer iruse(maxruns)
       common/lcdat1/nobs,p,v,ve
C       common/lcdat2/nslc,sp,sv
       common/lcdat2/nslc,ilc,sp,sv,iverb

C  PARAMETERS FOR ERROR TREATMENTS: FLOOR TO ERROR (WT).
       parameter (emin=0.005)

C  UNIT 9 HAS STD LCS
C       slcfile = 'stdlc.5tmplts'
C       slcfile = 'stdlc.3rr1sin2ecl'
       slcfile = 'stdlc.all10'
       write(6,90) slcfile
 90    format('Reading STD LCs from: ',a30)
       open(9,file=slcfile,status='old')
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  READ IN FILE CONTAINING NUMBER AND DATES OF OBSERVING RUNS (UNIT 17)
       open(17,file='stdlc_per10.runs', status='old')
       read(17,*) nruns
       write(6,672) nruns
 672   format('Number of observing runs = ',i2)
       do 673 i=1,nruns
          read(17,*) tbeg(i), tend(i)
          write(6,676) i, tbeg(i), tend(i)
 676      format('Run(i), Time(begin), Time(end) = ',i2,2f8.1)
          write(6,678)
 678      format('   Use data from this run? (1=yes, 0=no):  ',$)
          read(5,*) iruse(i)
 673      continue
C       write(6,921)
C 921   format('Observing runs: 94jun 94jun2 94oct 95may 96jun')
C       write(6,922)
C 922   format('Pick obsruns to use (1=use, 0=dont), eg  1 1 1 0 0:  ',$)
C       read(5,*) ir1,ir2,ir3,ir4,ir5

          
C  UNIT 8 HAS INPUT DATA FILE: READ IT
       write(6,92)
 92    format('Enter name of file w/ HJD Vmag Imag (vfe.B01):  ',$)
       read(5,100) infile
       write(6,923)
 923   format('Enter error ratio for bad seeing points (eg. 2.0):  ',$)
       read(5,*) erat
       write(6,924)
 924   format('Enter maximum error value to accept (eg. 0.1):  ',$)
       read(5,*) errmax
       open(8,file=infile,status='old')
       nobs = 0
       vmin = 99.
       do 10 i=1,99999
          read(8,*,err=10,end=11) time,iptp,vm,tverr,bm,be
          if(tverr.gt.1.0) goto 10
          if(tverr.gt.errmax) goto 10
C   DETERMINE WHETHER DATA SHOULD(NOT) BE USED -- IF IT FALLS W/IN
C     TBEG ADN TEND FOR A CERTAIN RUN, AND IRUSE=0 FOR THAT RUN,
C     DO NOT USE THAT DATA.
          do 483 k=1,nruns
             if(time.gt.tbeg(k) .and. time.lt.tend(k)) then
               if(iruse(k) .eq. 0) goto 10
             endif
 483         continue
CC   SELECT HJDS... HJDS OF DIFFERENT RUNS ARE HARDWIRED HERE...
CC     94JUN1
C          if(time.gt.520.2.and.time.lt.522.2) then
C            if(ir1.eq.0) goto 10
C          endif
CC     94JUN2
C          if(time.gt.529.2.and.time.lt.531.2) then
C            if(ir2.eq.0) goto 10
C          endif
CC     94OCT
C          if(time.gt.638.2.and.time.lt.642.2) then
C            if(ir3.eq.0) goto 10
C          endif
CC     95MAY
C          if(time.gt.896.2.and.time.lt.898.2) then
C            if(ir4.eq.0) goto 10
C          endif
CC     96JUN
C          if(time.gt.1278.2.and.time.lt.1279.2) then
C            if(ir5.eq.0) goto 10
C          endif
          nobs = nobs+1
          hjd(nobs) = time
          vmag(nobs) = vm
          imag(nobs) = bm
C    THIS PUTS A FLOOR ON THE ERROR, SO THE WEIGHTING DOESN'T GO NUTS
          if(tverr.lt.emin) tverr=emin
C    THIS ALLOWS THE BAD-SEEING DATA TO HAVE LESS INFLUENCE, BAD SEEING
C     SM PTYPES ALL END IN 0 (OPEN SYMBOLS).
          open = 10.0*(float(iptp)/10.0 - float(iptp/10))
          if(abs(open).lt.0.01) tverr = tverr*erat
C          WRITE(6,*) IPTP,OPEN,TVERR
          verr(nobs) = tverr
          if(vm.lt.vmin) then
            vmin = vm
            t0 = time
          endif
 10       continue
 11    write(6,78) nobs
 78    format(/,'READ ',i3,' mags from variable star data file',/)

C  UNIT 20 HAS FILE WITH OUTPUT CHISQ VALUES FOR EACH FIT.
       write(6,96)
 96    format('Enter name for output file:  ',$)
       read(5,100) outfile
       open(20,file=outfile,status='new')

100    format(a)

       write(6,64)
 64    format('Output verbosity in reporting indiv solns (hi/lo/no = 2/1/0):  ',$)
       read(5,*) iverb

C  LOOP 1: STEP THROUGH SERIES OF PERIODS AND DO FITS FOR EACH PER
       write(6,70)
 70    format('Enter Starting and Ending Periods, and Fractional Increment (eg 0.2 1.0 0.001):')
       read(5,*) per1,per2,dper
       chimax = -99.9
       chisum = 0.
       nchi = 0
       perold = per1
       do 20 i=1,99999
          per = (1.0+dper)*perold
          if(per.gt.per2) goto 222
          bobo = float(i)/300.0 - float(i/300)
          if(abs(bobo).lt.0.0001) write(6,77) i,per
 77       format('i  per(i) = ',i5,f11.7)
C  GIVEN THIS PERIOD, COMPUTE PHASES FOR EACH DATA POINT
          do 15 j=1,nobs
             cyc = (hjd(j)-t0)/per
             p(j) = cyc - float(int(cyc))
             if(p(j).lt.0.) p(j) = p(j)+1.0
             if(p(j).ge.1.) p(j) = p(j)-1.0
C             WRITE(6,79) J,HJD(J),CYC,P(J),VMAG(J),PER
C 79          FORMAT('j hjd(j) cyc p(j) = ',i3,2f12.4,3f8.4)
 15          continue
C  NOW FIT THE VBAND LIGHT CURVE WITH EACH TEMPLATE
C    CURRENTLY HARDWIRED FOR 5 V-ONLY TEMPLATES, BUT EXPANDABLE
          do 16 j=1,nobs
             v(j) = vmag(j)
             ve(j) = verr(j)
 16          continue
          ilc = 1
          call fitlc(xx,xsq)
C     THIS CHECKS FOR NaN ERRORS, IF PRESENT, SET TO 78.9
          if(xsq.ne.xsq) xsq = 78.9
          xsq1v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms1v = sqrt(xsq)
          ilc = 2
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq2v = xsq
          if(xsq.lt.10.0) then 
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms2v = sqrt(xsq)
          ilc = 3
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq3v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms3v = sqrt(xsq)
          ilc = 4
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq4v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms4v = sqrt(xsq)
          ilc = 5
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq5v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms5v = sqrt(xsq)
          ilc = 6
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq6v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms6v = sqrt(xsq)
          ilc = 7
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq7v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms7v = sqrt(xsq)
          ilc = 8
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq8v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms8v = sqrt(xsq)
          ilc = 9
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq9v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms9v = sqrt(xsq)
          ilc = 10
          call fitlc(xx,xsq)
          if(xsq.ne.xsq) xsq = 78.9
          xsq0v = xsq
          if(xsq.lt.10.0) then
            chimax = max(chimax,xsq)
            chisum = chisum+xsq
            nchi = nchi+1
          endif
C          rms0v = sqrt(xsq)
C KILL THIS PART IN THIS VERSION
CC  NOW FIT THE IBAND LIGHT CURVE WITH EACH TEMPLATE
C          do 17 j=1,nobs
C             v(j) = imag(j)
C 17          continue
C          ilc = 1
C          call fitlc(xx,xsq)
C          xsq1i = xsq
CC          rms1i = sqrt(xsq)
C          ilc = 2
C          call fitlc(xx,xsq) 
C          xsq2i = xsq
CC         rms2i = sqrt(xsq)
C          ilc = 3
C          call fitlc(xx,xsq)
C          xsq3i = xsq
CC          rms3i = sqrt(xsq)
C  WRITE THE RMS VALUES TO THE OUTPUT FILE AND SET PEROLD
          write(20,234) per,xsq1v,xsq2v,xsq3v,xsq4v,xsq5v,xsq6v,xsq7v,xsq8v,xsq9v,xsq0v
C          write(20,234) per,rms1v,rms2v,rms3v,rms1i,rms2i,rms3i
 234      format(f9.7,10f10.6)
          perold = per
          nperpts = nperpts+1
 20       continue

C   READ LAST FDF ... DONE WITH THIS SIMULATION
 222   write(6,2129) outfile
 2129  format(/,'Ending simulation, see ',a40)
       write(6,2229)
 2229  format(/,'Output = Period xsq(1) xsq(2) xsq(3) xsq(4) xsq(5) xsq(6) xsq(7) xsq(8) xsq(9) xsq(10)')
C 2229  format(/,'Output = Period xsq(1,V) xsq(2,V) xsq(3,V) xsq(1,I) xsq(2,I) xsq(3,I)')
C 2229  format(/,'Output = Period rms(1,V) rms(2,V) rms(3,V) rms(1,I) rms(2,I) rms(3,I)')
       write(6,3468) nperpts
 3468  format('N Per Pts = ',i5)
C 3468  format('N Per Pts = ',i6)
       chiave = chisum/float(nchi)
       chilim = 1.4*chiave
       write(6,3478) chimax,chilim
 3478  format(/,'Maximum chi^2 and suggested plot limit = ',2f11.4)

C  NOW MAKE SM PLOT MACRO TO PLOT THIS MONKEY -- HARDWIRED FOR 6 PLOTS,
C     BUT THATS OK SINCE FIRST 6 TMPLTS ARE THE INDIV RRABS, AND WE'VE
C     ALREADY LOOKED AT THE RRC,SIN,WUMA,ALGOL CURVES.
       open(21,file='sm.'//outfile,status='new')
       call mkpltsm(nperpts,per1,per2,chilim,outfile)
       write(6,4523)
 4523  format(/'To plot with SM:')
       write(6,4525) outfile
 4525  format(' macro read sm.',a20)
       write(6,4527)
 4527  format(' plt',/)

       close(21)
       close(20)
       close(8)
       close(9)
       stop
       end



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

      real p(256),v(256),sp(256),sv(10,256),ve(256)
      real pp(4,3),yy(4),xx(3)
      character*1 ckey
       common/lcdat1/nobs,p,v,ve
C       common/lcdat2/nslc,sp,sv
       common/lcdat2/nslc,ilc,sp,sv,iverb
 
C  IVERB=1 GIVES REAL TIME SCREEN OUTPUT OF FITTING, IVERB=0 DOESNT.
C      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
C         WRITE(6,*) I,I,P(I),V(I)
         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.2) 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)
C         WRITE(6,*) I,J,YY(I)
 32      continue

C  REALITY CHECK
      if(iverb.eq.2) 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.2) write(6,320) 
 320  format('          N       a        p0       Vi')
      do 40 i=1,4
         if(iverb.eq.2) 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.2) 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*BEST
      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.2) write(6,320) 
      do 70 i=1,4
         if(iverb.eq.2) 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  KLUDGE: IF ITER > ITMAX IN AMOEBA (IE, IT DIDN'T CONVERGE), SET 
C   XSQ LARGE SO IT DOESN'T CHOSE THIS FIT!
      if(iter.gt.4999) xsq = 99.99

C  FINAL VECTOR (THIS *HAS* BEEN RESTARTED NEAR FINAL SOLUTION)
      if(iverb.eq.1) 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)

C  KLUDGE2: IF XSQ < 1e-8, ITS A DEGENERATE FIT, SO DON'T USE IT.
      if(xsq.lt.1.0e-8) then
        write(6,719)
 719    format('WARNING: x^2 tiny (degenerate fit), so setting x^2=88.88')
        xsq = 88.88
      endif

      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)
      common/lcdat1/nobs,p,v,ve
C      common/lcdat2/nslc,sp,sv
      common/lcdat2/nslc,ilc,sp,sv,iverb

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)*ve(i))
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
      PARAMETER (NMAX=20,ITMAX=999)
C      PARAMETER (NMAX=20,ITMAX=50000)
      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     SUBROUTINE TO AUTOMATICALLY MAKE SM PLOT FILES
C       NB: CURRENTLY HARDWIRED FR 6 PLOTS:
C            1 = TEMPLATE2 (RRA2)
C            2 = TEMPLATE5 (RRB2)
C            3 = TEMPLATE7 (RRC)
C            4 = TEMPLATE8 (SINE)
C            5 = TEMPLATE9 (WUMA)
C            6 = TEMPLATE10(ALGOL)

       subroutine mkpltsm(nperpts,perlo,perhi,chilim,outfile)

       character*50 outfile

       write(21,200) 
 200   format('plt    erase')
       write(21,210) outfile 
 210   format('       data ',a20)
       write(21,220) nperpts
 220   format('       lines 1 ',i6)
       write(21,230) 
C 230   format('       read {per 1  xv1 2  xv2 3  xv3 4  xv4 5  xv5 6  xv6 7}')
 230   format('       read {per 1  xv1 3  xv2 6  xv3 8  xv4 9  xv5 10  xv6 11}')
       write(21,240) 
 240   format('       expand 1.01')
       write(21,250) 
 250   format('       relocate (7000 31500)')
       write(21,260) outfile
 260   format('       putlabel 6 ',a20)

C  NOW DO TOP OF 6 PLOTS
       write(21,1000) 
 1000  format('       location 4000 31000 26500 31000')
       write(21,1010) perlo, perhi, chilim
 1010  format('       limits ',2f5.1,' 0. ',f4.2)
C       write(21,1010) chilim
C 1010  format('       limits 0.2 1.0 0. ',f4.2)
       write(21,1020) chilim/2.0
 1020  format('       ticksize 0.02 0.1 0.01 ',f4.2)
       write(21,1030) 
 1030  format('       box 0 1')
       write(21,1040) 
 1040  format('       connect per xv1')
       write(21,1050) 
 1050  format('       relocate  (1500 28750)')
       write(21,1060) 
 1060  format('       putlabel 5 chi^2(T2)')
       write(21,1052) 
 1052  format('       relocate  (1500 27750)')
       write(21,1062) 
 1062  format('       putlabel 5 (RRa_2)')

C  NOW DO 2nd OF 6 PLOTS
       write(21,2000) 
 2000  format('       location 4000 31000 22000 26500')
       write(21,2010) perlo, perhi, chilim
 2010  format('       limits ',2f5.1,' 0. ',f4.2)
C       write(21,2010) chilim
C 2010  format('       limits 0.2 1.0 0. ',f4.2)
       write(21,2020) chilim/2.0
 2020  format('       ticksize 0.02 0.1 0.01 ',f4.2)
       write(21,2030) 
 2030  format('       box 0 1')
       write(21,2040) 
 2040  format('       connect per xv2')
       write(21,2050) 
 2050  format('       relocate  (1500 24250)')
       write(21,2060) 
 2060  format('       putlabel 5 chi^2(T5)')
       write(21,2052) 
 2052  format('       relocate  (1500 23250)')
       write(21,2062) 
 2062  format('       putlabel 5 (RRb_2)')

C  NOW DO 3rd OF 6 PLOTS
       write(21,3000) 
 3000  format('       location 4000 31000 17500 22000')
       write(21,3010) perlo, perhi, chilim
 3010  format('       limits ',2f5.1,' 0. ',f4.2)
C       write(21,3010) chilim
C 3010  format('       limits 0.2 1.0 0. ',f4.2)
       write(21,3020) chilim/2.0
 3020  format('       ticksize 0.02 0.1 0.01 ',f4.2)
       write(21,3030) 
 3030  format('       box 0 1')
       write(21,3040) 
 3040  format('       connect per xv3')
       write(21,3050) 
 3050  format('       relocate  (1500 19750)')
       write(21,3060) 
 3060  format('       putlabel 5 chi^2(T7)')
       write(21,3052) 
 3052  format('       relocate  (1500 18750)')
       write(21,3062) 
 3062  format('       putlabel 5 (RRc)')

C  NOW DO 4th OF 6 PLOTS
       write(21,4000) 
 4000  format('       location 4000 31000 13000 17500')
       write(21,4010) perlo, perhi, chilim
 4010  format('       limits ',2f5.1,' 0. ',f4.2)
C       write(21,4010) chilim
C 4010  format('       limits 0.2 1.0 0. ',f4.2)
       write(21,4020) chilim/2.0
 4020  format('       ticksize 0.02 0.1 0.01 ',f4.2)
       write(21,4030) 
 4030  format('       box 0 1')
       write(21,4040) 
 4040  format('       connect per xv4')
       write(21,4050) 
 4050  format('       relocate  (1500 15250)')
       write(21,4060) 
 4060  format('       putlabel 5 chi^2(T8)')
       write(21,4052) 
 4052  format('       relocate  (1500 14250)')
       write(21,4062) 
 4062  format('       putlabel 5 (Sine)')

C  NOW DO 5th OF 6 PLOTS
       write(21,5000) 
 5000  format('       location 4000 31000 8500 13000')
       write(21,5010) perlo, perhi, chilim
 5010  format('       limits ',2f5.1,' 0. ',f4.2)
C       write(21,5010) chilim
C 5010  format('       limits 0.2 1.0 0. ',f4.2)
       write(21,5020) chilim/2.0
 5020  format('       ticksize 0.02 0.1 0.01 ',f4.2)
       write(21,5030) 
 5030  format('       box 0 1')
       write(21,5040) 
 5040  format('       connect per xv5')
       write(21,5050) 
 5050  format('       relocate  (1500 10750)')
       write(21,5060) 
 5060  format('       putlabel 5 chi^2(T9)')
       write(21,5052) 
 5052  format('       relocate  (1500 9750)')
       write(21,5062) 
 5062  format('       putlabel 5 (W-UMa)')

C  NOW DO BOTTOM OF 6 PLOTS
       write(21,6000) 
 6000  format('       location 4000 31000 4000 8500')
       write(21,6010) perlo, perhi, chilim
 6010  format('       limits ',2f5.1,' 0. ',f4.2)
C       write(21,6010) chilim
C 6010  format('       limits 0.2 1.0 0. ',f4.2)
       write(21,6020) chilim/2.0
 6020  format('       ticksize 0.02 0.1 0.01 ',f4.2)
       write(21,6030) 
 6030  format('       box 1 1')
       write(21,6040) 
 6040  format('       connect per xv6')
       write(21,6050) 
 6050  format('       relocate  (1500 6250)')
       write(21,6060) 
 6060  format('       putlabel 5 chi^2(T10)')
       write(21,6052) 
 6052  format('       relocate  (1500 5250)')
       write(21,6062) 
 6062  format('       putlabel 5 (Algol)')
       write(21,6070) 
 6070  format('       xlabel Period [days]')

       return
       end
