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 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 , 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 o Vmaxo Vmino Amplo dPro Emaxo Nobs 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 o Vmaxo Vmino Amplo dPro Emaxo Nobs 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 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