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