# File rvsao/Xcsao/xcfit.x # March 11, 1999 # By Doug Mink, Harvard-Smithsonian Center for Astrophysics # After Gerard Kriss, Johns Hopkins University and others # Copyright(c) 1999 Smithsonian Astrophysical Observatory # You may do anything you like with this file except remove this copyright. # The Smithsonian Astrophysical Observatory makes no representations about # the suitability of this software for any purpose. It is provided "as is" # without express or implied warranty. # XCFIT finds the best template and velocity for one object spectrum include include include include include "rvsao.h" include "xcv.h" procedure xcfit (specfile,specdir,mspec,mband) char specfile[ARB] # Object spectrum file name char specdir[ARB] # Directory for object spectra int mspec # Object aperture to read from multispec file int mband # Object band to read from multispec file pointer spectrum # Object spectrum pointer specim # Object image header structure bool arcflag # archive record writing flag bool newfit # true to redo cross-correlation int i, j char specpath[SZ_PATHNAME] # Object spectrum path name char tempfiles[SZ_PATHNAME] # List of template spectra char tempfile[SZ_PATHNAME] # Template spectrum file name char temppath[SZ_PATHNAME] # Template spectrum path name char tempdir[SZ_PATHNAME] # Directory for template spectra char tempnums0[SZ_LINE] # List of multispec spectra to read char tempnums[SZ_LINE] # List of multispec spectra to read char txfile[SZ_PATHNAME,MAXTEMPS] # Correlation file names int rmode # Report format (1=normal,2=one-line) double tz # Corrections to correlation Z char xlab[SZ_LINE] # Title for wavelength plots of spectrum char chopstr[SZ_LINE] # Chop out emission lines # (yes | no | tempfile | specfile) char zero_pad[SZ_LINE] # Zero-pad before correlating or from file # (no | yes | tempfile) pointer specwcs # Pointer to first data point with WCS pointer wlwcs # Pointer to first wavelength with WCS int npspec # Number of spectrum points with WCS pointer tempwcs # Pointer to first template data point with WCS pointer twlwcs # Pointer to first template wavelength with WCS int nptemp # Number of template points with WCS pointer xcor # cross-correlation returned from xcorfit pointer ixcor # cross-correlation returned from xcorfit pointer xvel # cross-correlation velocities from xcorfit pointer ixvel # cross-correlation velocities from xcorfit int ncvel[MAXTEMPS] # Number of points in each correlation double zshift # guessed Z for template offset double czguess # initial velocity or Z for template offset double tmstrt,tmfnsh # template limits after zguess applied bool renorm # renormalization flag int izpass # Count for template shifts int nzpass # Number of template shifts bool zinit # TRUE if initial value for z is set int tspec # Template spectrum to read from multispec file int tzpad # Flag for zero-padding for one template int ntcor # ntemp * ncor2 bool chop0 # TRUE if emission lines already removed from spectrum pointer templist # List of template files double tcz # Velocity returned from xcorfit double tczerr # Velocity error returned from xcorfit double tczr # R returned from xcorfit double rmax # Maximum R value for spectrum and best template pointer shspec # Spectrum pixels in log-lambda for xcorfit pointer shtemp # Template pixels in log-lambda for xcorfit pointer wltemp # Wavelength vector for template overlap double minwav0,maxwav0 # Starting,ending wavelengths from parameter file double minwav1,maxwav1 # Overlap wavelength limits from spectrum and template double minwav,maxwav # Wavelength limits for cross-correlation double logmin,logmax # Log-10 of wavelength limits for cross-correlation bool filexcor # Flag to save correlation to a file double spwl1,spwl2 # Spectrum first and last wavelengths in angstroms double tmwl1,tmwl2 # Template first and last wavelengths in angstroms int ntx # Number of templates actually cross-correlated int tband # Multispec band for templates real spmean, speci double sptot double spshft, tmshft # Spectrum shift in log wavelength double maxstrt, minfin, dlog10() real emrej1,emrej2,rejem,abrej1,abrej2,rejabs int itemp, nsp, ncor, ncor2 double tmpvel pointer tempspec # Template spectrum pointer tempim # Template image header structure pointer tempsh # Template spectrum header structure int ntspec,ntspec0 # Number of template multispec spectra int tspec_range[3,MAX_RANGES] bool echelle # If true, template multispec numbers track object int ip,jp,lfile # Limits for multispec aperture decoding char lbracket[3] # left bracket, parenthesis, brace char rbracket[3] # right bracket, parenthesis, brace bool correlate # True to cross-correlate spectrum against templates double dindef char cdot bool clgetb() int clgeti(), clscan() real clgetr(), max(), min() double clgetd() int strdic(), stridx(), stridxs() int decode_ranges(),get_next_number() int imtgetim(), imaccess(), strlen(), ldir pointer imtopen() pointer sp include "rvsao.com" include "results.com" include "contin.com" define refit_ 30 define newtemp_ 40 define newtap_ 50 define endspec_ 60 define endcorr_ 70 define endxc_ 90 begin call smark (sp) c0 = 299792.5 cdot = '.' dindef = INDEFD specsh = NULL tempsh = NULL call sprintf (lbracket,3,"[({") call sprintf (rbracket,3,"])}") emrej1 = clgetr ("s_em_reject") abrej1 = clgetr ("s_abs_reject") emrej2 = clgetr ("t_em_reject") abrej2 = clgetr ("t_abs_reject") correlate = TRUE correlate = clgetb ("correlate") if (!correlate) savevel = FALSE newresults = FALSE newfit = FALSE # Load spectrum ldir = strlen (specdir) if (ldir > 0) { if (specdir[ldir] != '/') { specdir[ldir+1] = '/' specdir[ldir+2] = EOS } call strcpy (specdir,specpath,SZ_PATHNAME) call strcat (specfile,specpath,SZ_PATHNAME) } else call strcpy (specfile,specpath,SZ_PATHNAME) # Check for readability of object spectrum if (imaccess (specpath, READ_ONLY) == NO) { call eprintf ("XCFIT: cannot read spectrum path %s \n") call pargstr (specpath) call sfree (sp) return } call getspec (specpath, mspec, mband, spectrum, specim, savevel) if (specim == ERR) { # call eprintf ("XCFIT: Error reading spectrum %s\n") # call pargstr (specpath) go to endxc_ } call strcpy (LABEL(specsh), xlab, SZ_LINE) if (strlen (xlab) <= 0) call strcpy ("Wavelength", xlab, SZ_LINE) if (strlen (UNITS(specsh)) > 0) { call strcat (" in ",xlab,SZ_LINE) call strcat (UNITS(specsh),xlab,SZ_LINE) } else call strcat (" in Angstroms",xlab,SZ_LINE) # Renormalize spectrum if requested renorm = clgetb ("renormalize") if (renorm) { sptot = 0.d0 nsp = 0 do i = 1, specpix { speci = Memr[spectrum+i-1] if (speci != 0.) { sptot = sptot + speci nsp = nsp + 1 } } if (nsp > 0) { spmean = 0.001d0 * sptot / double (nsp) do i = 1, specpix { if (Memr[spectrum+i-1] != 0.) { Memr[spectrum+i-1] = Memr[spectrum+i-1] / spmean } } } else { call eprintf ("XCFIT: Spectrum is all zeroes\n") go to endxc_ } } # Set pixel limits of spectrum WCS specwcs = spectrum + NP1(specsh) - 1 wlwcs = SX(specsh) + NP1(specsh) - 1 npspec = NP2(specsh) - NP1(specsh) + 1 # If plot enabled, show the object spectrum. if (pltspec) call plotspec (npspec, Memr[specwcs], specname, Memr[wlwcs], xlab, nsmooth) # Eliminate bad lines from spectrum if (clscan("fixbad") != EOF) { if (clgetb ("fixbad")) { call filllist (npspec, Memr[specwcs], Memr[wlwcs],debug) if (pltspec) call plotspec (npspec,Memr[specwcs], specname,Memr[wlwcs],xlab,nsmooth) } } # Optional wavelength limits for cross-correlation minwav0 = clgetd ("st_lambda") if (minwav0 == 0.d0) minwav0 = dindef maxwav0 = clgetd ("end_lambda") if (maxwav0 == 0.d0) maxwav0 = dindef # Templates against which to correlate spectra call clgstr ("templates",tempfiles,SZ_PATHNAME) call clgstr ("tempdir",tempdir,SZ_PATHNAME) ldir = strlen (tempdir) if (ldir > 0 && tempdir[ldir] != '/' && tempdir[ldir] != '$') { tempdir[ldir+1] = '/' tempdir[ldir+2] = EOS } # Multispec template numbers (use only first if multiple files) tband = clgeti ("tempband") call clgstr ("tempnum",tempnums0,SZ_LINE) if (decode_ranges (tempnums0, tspec_range, MAX_RANGES, ntspec0) == ERR) call error (1, "XCFIT: Illegal multispec list") # Echelle switch (template spectrum tracks object spectrum number if true) echelle = clgetb ("echelle") # Zero padding on, off, or read from template file call clgstr ("zeropad",zero_pad,SZ_LINE) tzpad = strdic (zero_pad,zero_pad,SZ_LINE, XC_ZTYPES) # Emission line chopping flag call clgstr ("s_emchop",chopstr,SZ_LINE) sfchop = strdic (chopstr, chopstr, SZ_LINE, XC_CTYPES) call clgstr ("t_emchop",chopstr,SZ_LINE) tfchop = strdic (chopstr, chopstr, SZ_LINE, XC_CTYPES) # Fraction of spectrum to apodize at each end han = clgetr ("bell_window") # Filter type and filter constants [always cosine-bell] filter = COSBELL lo = clgeti ("low_bin") toplo = clgeti ("top_low") nrun = clgeti ("nrun") topnrn = clgeti ("top_nrun") if (debug) { call printf ("XCFIT: filter is %d %d %d %d\n") call pargi (lo) call pargi (toplo) call pargi (topnrn) call pargi (nrun) } # Optional tshift to be used for all spectra and templates # tshift can also be read from the TSHIFT parameter in a template header tshift = clgetd ("tshift") # If npts not a power of two, determine next largest power of two: npts = clgeti ("ncols") if (npts != 256 && npts != 512 && npts != 1024 && npts != 2048 && npts != 4096 && npts != 8192 && npts != 16384) { npts = alog (real (npts)) / LN2 + 1. npts = 2 ** npts } # Check for unreasonable filtering if (lo < 1) lo = 1 if (lo > npts) lo = 1 if (toplo < 1) toplo = 1 if (toplo > npts) toplo = 1 if (nrun < 1) nrun = npts if (nrun > npts) nrun = npts if (topnrn < 1) topnrn = npts if (topnrn > npts) topnrn = npts # Allocate correlation and velocity vectors ncor2 = npts * 2 call countemp (tempfiles,tempdir,tempnums0,mspec,echelle,debug,itemp) if (debug) { call printf ("XCFIT: %d points in spectrum, up to %d points in transforms\n") call pargi (npts) call pargi (ncor2) call printf ("XCFIT: %d templates being cross-correlated\n") call pargi (itemp) call flush (STDOUT) } ntcor = ncor2 * itemp call salloc (xcor, ntcor, TY_REAL) call salloc (xvel, ntcor, TY_REAL) # Optional object velociity to be used to shift templates czguess = 0.d0 zinit = FALSE switch (vinit) { case IVZERO: czguess = 0.d0 zinit = FALSE case IVGUESS: czguess = clgetd ("czguess") zinit = TRUE case IZGUESS: czguess = c0 * clgetd ("czguess") zinit = TRUE case IVXC: if (spxvel != dindef) { czguess = spxvel zinit = TRUE } case IVEM: if (spevel != dindef) { czguess = spevel zinit = TRUE } case IVCOMB: if (spvel != dindef) { czguess = spvel zinit = TRUE } default: } # Number of times to shift template nzpass = 0 nzpass = clgeti ("nzpass") nzpass = nzpass + 1 # If npts not a power of two, determine next largest power of two: npts = clgeti ("ncols") if (npts != 256 && npts != 512 && npts != 1024 && npts != 2048 && npts != 4096 && npts != 8192 && npts != 16384) { npts = alog (real (npts)) / LN2 + 1. npts = 2 ** npts } # Allocate rebinned spectra and wavelength vectors call salloc (shspec, npts, TY_REAL) call salloc (shtemp, npts, TY_REAL) call salloc (wltemp, npts, TY_REAL) # If not cross-correlating spectrum skip to result plotting if (!correlate) { if (ntemp < 1) ntemp = 1 ntx = ntemp if (itmax < 1) itmax = 1 minwav = W0(specsh) if (minwav0 != dindef && minwav0 > minwav) minwav = minwav0 maxwav = W1(specsh) if (maxwav0 != dindef && maxwav0 < maxwav) maxwav = maxwav0 do itemp = 1, ntemp { twl1[itemp] = minwav twl2[itemp] = maxwav } go to endcorr_ } # For each template, compute cross-correlation refit_ templist = imtopen (tempfiles) rmax = 0.d0 ntx = 0 itemp = 0 ntemp = 0 spvqual = 0 # Initialize new template file newtemp_ tspec = -1 # Get next template spectrum file name from the list if (imtgetim (templist, tempfile, SZ_PATHNAME) == EOF) go to endspec_ # Check for specified apertures in multispec template file ip = stridxs (lbracket,tempfile) if (echelle) ntspec = 1 else if (ip > 0) { lfile = strlen (tempfile) tempfile[ip] = EOS ip = ip + 1 jp = 1 while (stridx (tempfile[ip],rbracket) == 0 && ip <= lfile) { tempnums[jp] = tempfile[ip] tempfile[ip] = EOS ip = ip + 1 jp = jp + 1 } tempnums[jp] = EOS if (decode_ranges (tempnums,tspec_range,MAX_RANGES,ntspec) == ERR) call error (1, "XCFIT: Illegal template multispec list") } else { call strcpy (tempnums0,tempnums,SZ_LINE) if (decode_ranges (tempnums,tspec_range,MAX_RANGES,ntspec) == ERR) call error (1, "XCFIT: Illegal template multispec list") } # if (debug) { # call printf ("XCFIT: next template file is %s [%s] = %d aps\n") # call pargstr (tempfile) # call pargstr (tempnums) # call pargi (ntspec) # } # Check for readability of template spectrum if (stridx ("/",tempfile) == 1 || stridx ("$",tempfile) > 0 || stridx (cdot,tempfile) == 1) call strcpy (tempfile,temppath,SZ_PATHNAME) else { call strcpy (tempdir,temppath,SZ_PATHNAME) call strcat (tempfile,temppath,SZ_PATHNAME) } if (imaccess (temppath, READ_ONLY) == NO) { call eprintf ("XCFIT: cannot read template file %s \n") call pargstr (temppath) go to newtemp_ } # Get next multispec aperture number from list newtap_ if (echelle) tspec = mspec else if (get_next_number (tspec_range, tspec) == EOF) go to newtemp_ # Initialize cross-correlation results for this template itemp = itemp + 1 tz = 0.d0 tempvel[itemp] = 0.d0 tempshift[itemp] = 0.d0 temphcv[itemp] = 0.d0 cz[itemp] = 0.d0 czerr[itemp] = 0.d0 czr[itemp] = 0.d0 zvel[itemp] = 0.d0 # Load template spectrum call gettemp (tempfile,tspec,tband,tempdir,tempspec,tempim,tempsh,itemp) if (tempim == ERR) go to newtap_ ntx = ntx + 1 ntemp = itemp if (tspec > 0) { call sprintf (tempid[1,itemp], SZ_PATHNAME,"%s[%d]") call pargstr (tempfile) call pargi (tspec) } else { call sprintf (tempid[1,itemp], SZ_PATHNAME, "%s") call pargstr (tempfile) } # Set pixel limits of template WCS tempwcs = tempspec + NP1(tempsh) - 1 twlwcs = SX(tempsh) + NP1(tempsh) - 1 nptemp = NP2(tempsh) - NP1(tempsh) + 1 if (plttemp) call plotspec (nptemp,Memr[tempwcs],tempfile, Memr[twlwcs],xlab,nsmooth) # Spectrum emission line chopping flag schop = FALSE emrej[1] = emrej1 abrej[1] = abrej1 switch (sfchop) { case CHOP: schop = TRUE case NOCHOP: schop = FALSE case TEMPFILE: schop = TRUE call imgbpar (tempim,"CHOPEM ",schop) # Read spectrum chopping limits from template header, if present if (schop) { rejem = emrej[1] call imgrpar (tempim,"EMREJ",rejem) emrej[1] = rejem rejabs = abrej[1] call imgrpar (tempim,"ABSREJ",rejabs) abrej[1] = rejabs if (debug) { call printf ("XCFIT: Keep spectrum data from -%f to %f sigma\n") call pargr (abrej[1]) call pargr (emrej[1]) } } case SPECFILE: schop = TRUE call imgbpar (specim,"CHOPEM ",schop) } chop0 = FALSE call imgbpar (specim,"EMCHOP ",chop0) if (chop0) schop = FALSE if (abrej[1] < emrej[1]) { tachop[itemp] = schop tschop[itemp] = FALSE } else { tschop[itemp] = schop tachop[itemp] = FALSE } # Template emission line chopping flag tchop = FALSE abrej[2] = abrej2 emrej[2] = emrej2 switch (tfchop) { case CHOP: tchop = TRUE case NOCHOP: tchop = FALSE case TEMPFILE: tchop = TRUE call imgbpar (tempim,"CHOPEM",tchop) # Read template chopping limits from template header, if present if (tchop) { rejem = emrej[2] call imgrpar (tempim,"EMREJ",rejem) emrej[2] = rejem rejabs = abrej[2] call imgrpar (tempim,"ABSREJ",rejabs) abrej[2] = rejabs if (debug) { call printf ("XCFIT: Keep template data from -%f to %f sigma\n") call pargr (abrej[2]) call pargr (emrej[2]) } } case SPECFILE: tchop = TRUE call imgbpar (specim,"CHOPEM",tchop) } chop0 = FALSE call imgbpar (tempim,"EMCHOP",chop0) if (chop0) tchop = FALSE if (debug) { if (tchop) call printf ("XCFIT: Template emission lines removed\n") else call printf ("XCFIT: Template emission lines left in\n") } # Turn off template continuum subtraction, if not needed tcont = TRUE call imgbpar (tempim,"SUBCONT",tcont) if (debug) { if (tcont) call printf ("XCFIT: Template continuum subtracted\n") else call printf ("XCFIT: Template continuum already subtracted\n") } tscont[itemp] = tcont # Divide spectrum and template by their continuum divcon = FALSE call imgbpar (tempim,"DIVCONT",divcon) if (debug) { if (divcon) call printf ("XCFIT: Continuum divided, not subtracted\n") else call printf ("XCFIT: Continuum subtracted, not divided\n") } tdivcon[itemp] = divcon # Ignore wavelength limit parameters and use file limits overlap = FALSE call imgbpar (tempim,"OVERLAP",overlap) if (minwav0 == dindef && maxwav0 == dindef) overlap = FALSE if (debug) { if (overlap) call printf ("XCFIT: Use entire overlapping wavelength space\n") else call printf ("XCFIT: Use wavelength limit parameters, if set\n") } toverlap[itemp] = overlap # Zero-padding of transforms before cross-correlation zpad = FALSE switch (tzpad) { case ZPAD: zpad = TRUE case NOZPAD: zpad = FALSE case TEMPFILE: call imgbpar (tempim,"ZEROPAD",zpad) } # Set number of points in correlation depending on whether zero-padding if (zpad) { ncor = 2 * npts if (debug) call printf ("XCFIT: Correlation zero-padded\n") } else { ncor = npts if (debug) call printf ("XCFIT: Correlation not zero-padded\n") } # Set fraction of correlation peak to fit pkfrac = clgetd ("pkfrac") call imgdpar (tempim,"PEAKFRAC",pkfrac) # Set redshift velocity limits minvel = clgetd ("minvel") maxvel = clgetd ("maxvel") if (maxvel != dindef && minvel != dindef && maxvel < minvel) { tmpvel = minvel minvel = maxvel maxvel = tmpvel } # Add known velocity components as 1+z's tz = (1.d0 + (spechcv / c0)) * (1.d0 + (tempvel[itemp] / c0)) * (1.d0 + (tempshift[itemp] / c0)) * (1.d0 + (tshift / c0)) / (1.d0 + (temphcv[itemp] / c0)) if (debug) { call printf ("XCFIT: spechcv= %.4f, tempvel= %.4f, temphcv= %.4f\n") call pargd (spechcv) call pargd (tempvel[itemp]) call pargd (temphcv[itemp]) call printf ("XCFIT: tz = %f\n") call pargd (tz) } # Set z offset for initial pass through loop if (zinit) { zshift = (1.d0 + (czguess / c0)) / tz if (debug) { call printf ("XCFIT: zshift is %f\n") call pargd (zshift) } } else zshift = 1.d0 # If shifting template on second pass, start here do izpass = 1, nzpass { tmstrt = testrt[itemp] tmfnsh = tefnsh[itemp] # Modify wavelength limits of template for velocity guess if (zinit || izpass > 1) { tmshft = dlog10 (zshift) tmstrt = tmstrt + tmshft tmfnsh = tmfnsh + tmshft } else tmshft = 0.d0 # Check to make sure that template and spectrum overlap if (spstrt > tmfnsh) { spwl1 = 10.d0 ** spstrt tmwl2 = 10.d0 ** tmfnsh call printf ("XCFIT: Spectrum starts at %9.2f angstroms\n") call pargd (spwl1) call printf ("XCFIT: Template ends at %9.2f angstroms\n") call pargd (tmwl2) go to endxc_ } if (spfnsh < tmstrt) { spwl2 = 10.d0 ** spfnsh tmwl1 = 10.d0 ** tmstrt call printf ("XCFIT: Spectrum ends at %9.2f angstroms\n") call pargd (spwl2) call printf ("XCFIT: Template starts at %9.2f angstroms\n") call pargd (tmwl1) go to endxc_ } # Figure lambda and log-lambda limits and increment maxstrt = max (spstrt, tmstrt) minwav1 = 10. ** maxstrt minfin = min (spfnsh, tmfnsh) maxwav1 = 10. ** minfin if (overlap) minwav = minwav1 else if ((minwav0==dindef)||(minwav0maxwav1)) maxwav = maxwav1 else maxwav = maxwav0 logmin = dlog10 (minwav) logmax = dlog10 (maxwav) twl1[itemp] = minwav twl2[itemp] = maxwav # Fill in globally available values of wavelength wave0 = minwav logw0 = logmin dlogw = (logmax - logmin) / double (npts - 1) deltav = dlogw * CLN10 if (debug) { call printf("XCFIT: overlap %9.3f to %9.3f by %9.3f, npts = %d \n") call pargd (10.d0 ** logmin) call pargd (10.d0 ** logmax) call pargd (10.d0 ** dlogw) call pargi (npts) } # Rebin spectrum to overlap this template if (debug) { call printf("XCFIT: spectrum %9.3f to %9.3f = %d %d\n") call pargd (10.d0 ** spstrt) call pargd (10.d0 ** spfnsh) call pargi (SN(specsh)) call pargi (DC(specsh)) } spshft = 0.d0 call rebinl (Memr[spectrum],specsh, spshft, Memr[shspec], npts, logw0, dlogw) # Rebin template to overlap this spectrum if (debug) { call printf("XCFIT: template %9.7f to %9.7f = %d %d\n") call pargd (10.d0 ** tmstrt) call pargd (10.d0 ** tmfnsh) call pargi (SN(tempsh)) call pargi (DC(tempsh)) } call rebinl (Memr[tempspec], tempsh, tmshft, Memr[shtemp], npts, logw0, dlogw) if (debug) { call printf ("XCFIT: Template %f %f %f\n") call pargr (Memr[tempspec]) call pargr (Memr[tempspec+1]) call pargr (Memr[tempspec+2]) call printf ("XCFIT: Template %f %f %f\n") call pargr (Memr[shtemp]) call pargr (Memr[shtemp+1]) call pargr (Memr[shtemp+2]) } # Set up wavelength vector for intermediate graphs do j = 1, npts { Memr[wltemp+j-1] = 10 ** (logw0 + (j-1) * dlogw) } # Set known portion of velocity if (zinit || izpass > 1) tvel = ((tz * zshift) - 1.d0) * c0 else tvel = (tz - 1.d0) * c0 if (debug) { call printf ("XCFIT: tvel = %f\n") call pargd (tvel) call printf ("XCFIT: filter is %d %d %d %d\n") call pargi (lo) call pargi (toplo) call pargi (topnrn) call pargi (nrun) } # Cross-correlate this template/spectrum combination ixcor = xcor + ((itemp - 1) * ncor2) ixvel = xvel + ((itemp - 1) * ncor2) call xcorfit (Memr[shspec],Memr[shtemp],Memr[wltemp], ncor, Memr[ixcor],Memr[ixvel],itemp, tcz,tczerr,tczr) # Save velocity and R-value cz[itemp] = tcz czerr[itemp] = tczerr czr[itemp] = tczr zvel[itemp] = (1.d0 + (tcz / c0)) * tz if (zinit || izpass > 1) { zvel[itemp] = zvel[itemp] * zshift tvshift[itemp] = (zshift - 1.d0) * c0 } else tvshift[itemp] = 0.d0 zvel[itemp] = (zvel[itemp] - 1.d0) * c0 if (nzpass > 1 && izpass < nzpass) { zshift = zshift * (1.d0 + (tcz / c0)) if (debug) { call printf ("XCFIT: pass %d: v= %9.3f, 1+z= %7.5f") call pargi (izpass) call pargd (tcz) call pargd (zshift) call printf (" r= %6.3f h= %6.4f arms= %8.6f\n") call pargd (tczr) call pargd (thght[itemp]) call pargd (tarms[itemp]) } } } # End of cz iteration if (debug) { call printf("End of cz iteration for template %d, %d points\n") call pargi (itemp) call pargi (ncor) call flush (STDOUT) } # Save velocity and cross-correlation vectors for summary plot ncvel[itemp] = ncor if (ncor < ncor2) { do i = ncor+1, ncor2 { Memr[ixcor+i] = 0.0 Memr[ixvel+i] = 0.0 } } if (debug) { call printf("Cross-correlation saved\n") call flush (STDOUT) } # if R is new max, set up velocity vector for correlation plot if (tczr >= rmax) { rmax = tczr npmax = tnpfit[itemp] xcrmax = zvel[itemp] itmax = itemp if (debug) { call printf("Best correlation template: %s %5.2f\n") call pargstr(tempname[1,itemp]) call pargd (rmax) } } newresults = TRUE # Print the filtered cross-correlation function if requested filexcor = FALSE filexcor = clgetb ("xcor_file") if (filexcor) call xcfile (specfile,specim,mspec,tempfile,tspec,itemp,ncor, Memr[ixvel],Memr[ixcor],txfile[1,itemp]) # Close the template image call close_image (tempim,tempsh) ntspec = ntspec - 1 if (ntspec > 0) go to newtap_ go to newtemp_ # Close the template list for this object spectrum endspec_ call imtclose (templist) endcorr_ # Save results for this object spectrum if at least one template correlated IM_UPDATE(specim) = NO if (ntx > 0) { spxvel = zvel[itmax] spxerr = czerr[itmax] spxr = czr[itmax] # Combine correlation and emission velocities for this object spectrum if (newresults) { call vcombine (spxvel,spxerr,spxr,spevel,speerr,spnlf, spvel,sperr,debug) # Save correlation information to archive file for this object arcflag = FALSE arcflag = clgetb ("archive") if (arcflag) call xcarch (specpath) } rmode = 2 rmode = clgeti ("report_mode") call xcrslts (specfile,mspec,specim,ncor,rmode,filexcor,txfile) # Plot results if (debug) { ixcor = xcor + ((itmax - 1) * ncor2) ixvel = xvel + ((itmax - 1) * ncor2) call printf ("XCFIT: template %d: vel= %.3f - %.3f\n") call pargi (itmax) call pargr (Memr[ixvel]) call pargr (Memr[ixvel+ncvel[itmax]-1]) call printf ("XCFIT: template %d: xcor= %.3f - %.3f\n") call pargi (itmax) call pargr (Memr[ixcor]) call pargr (Memr[ixcor+ncvel[itmax]-1]) } call xcplot (Memr[wlwcs],Memr[specwcs],npspec,Memr[xvel], Memr[xcor],ncvel,ncor2,specfile,mspec,specim,newfit) if (newfit) { newfit = FALSE go to refit_ } # Save best Cz, error, and R to image header, if requested # after checking for writability if (savevel || qplot) { if (imaccess (specpath, READ_WRITE) == NO) { call eprintf ("XCSAO: cannot write to %s; not saving results\n") call pargstr (specpath) } else if (savevel) { call vwhead (mspec,specim) call xcwhead (mspec,specim) } else if (qplot && newresults) { call qwhead (mspec,specim,"xcsao") call xcwhead (mspec,specim) } else if (qplot) { call qwhead (mspec,specim,"xcsao") } } } else { call printf ("XCFIT: No templates; no cross-correlation\n") } # Free rebinned spectra, wavelength, and correlation vectors endxc_ call sfree (sp) # Close the object spectrum image call close_image (specim, specsh) return end procedure countemp (tempfiles, tempdir, tempnums0, mspec, echelle, debug, ntemp) char tempfiles[ARB] # List of template spectra char tempdir[ARB] # Directory for template spectra char tempnums0[ARB] # List of multispec spectra to read int mspec # Object aperture to read from multispec file bool echelle # If true, template apertures track object bool debug int ntemp # Number of templates (returned) char tempfile[SZ_PATHNAME] # Template spectrum file name char temppath[SZ_PATHNAME] # Template spectrum path name char tempnums[SZ_LINE] # List of multispec spectra to read pointer templist # List of template files int ntspec,ntspec0 # Number of template multispec spectra int tspec_range[3,MAX_RANGES] int tspec # Template spectrum to read from multispec file int ip, jp, lfile char lbracket[3] # left bracket, parenthesis, brace char rbracket[3] # right bracket, parenthesis, brace char cdot int stridx(), stridxs() int decode_ranges(),get_next_number() int imtgetim(), imaccess(), strlen() pointer imtopen() define newtemp_ 40 define newtap_ 50 define endspec_ 60 begin call sprintf (lbracket,3,"[({") call sprintf (rbracket,3,"])}") cdot = '.' templist = imtopen (tempfiles) if (decode_ranges (tempnums0, tspec_range, MAX_RANGES, ntspec0) == ERR) call error (1, "COUNTEMP: Illegal multispec list") ntemp = 0 # Initialize new template file newtemp_ tspec = -1 # Get next template spectrum file name from the list if (imtgetim (templist, tempfile, SZ_PATHNAME) == EOF) go to endspec_ # Check for specified apertures in multispec template file ip = stridxs (lbracket,tempfile) if (echelle) ntspec = 1 else if (ip > 0) { lfile = strlen (tempfile) tempfile[ip] = EOS ip = ip + 1 jp = 1 while (stridx (tempfile[ip],rbracket) == 0 && ip <= lfile) { tempnums[jp] = tempfile[ip] tempfile[ip] = EOS ip = ip + 1 jp = jp + 1 } tempnums[jp] = EOS if (decode_ranges (tempnums,tspec_range,MAX_RANGES,ntspec) == ERR) call error (1, "COUNTEMP: Illegal template multispec list") } else { call strcpy (tempnums0,tempnums,SZ_LINE) if (decode_ranges (tempnums,tspec_range,MAX_RANGES,ntspec) == ERR) call error (1, "COUNTEMP: Illegal template multispec list") } if (debug) { call printf ("COUNTEMP: template file %d is %s [%s] = %d aps\n") call pargi (ntemp+1) call pargstr (tempfile) call pargstr (tempnums) call pargi (ntspec) call flush (STDOUT) } # Check for readability of template spectrum if (stridx ("/",tempfile) == 1 || stridx ("$",tempfile) > 0 || stridx (cdot,tempfile) == 1) call strcpy (tempfile,temppath,SZ_PATHNAME) else { call strcpy (tempdir,temppath,SZ_PATHNAME) call strcat (tempfile,temppath,SZ_PATHNAME) } if (imaccess (temppath, READ_ONLY) == NO) { call eprintf ("COUNTEMP: cannot read template file %s \n") call pargstr (temppath) go to newtemp_ } if (debug) { call printf ("COUNTEMP: template path %d is %s\n") call pargi (ntemp+1) call pargstr (temppath) call flush (STDOUT) } # Get next multispec aperture number from list newtap_ if (echelle) tspec = mspec else if (get_next_number (tspec_range, tspec) == EOF) go to newtemp_ # Add this template to counter ntemp = ntemp + 1 if (debug) { call printf ("COUNTEMP: %d templates\n") call pargi (ntemp) call flush (STDOUT) } ntspec = ntspec - 1 if (ntspec > 0) go to newtap_ go to newtemp_ endspec_ call imtclose (templist) if (debug) { call printf ("COUNTEMP: counted %d templates\n") call pargi (ntemp) call flush (STDOUT) } return end # Jul 13 1994 New subroutine # Jul 29 1994 Move header writing to xchead # Aug 3 1994 Fix erroneous subroutine calls # Aug 3 1994 Set emission line chopping switch for each spectrum template # Aug 8 1994 Add optional transform zero-padding # Aug 8 1994 Add spectrum directory # Aug 10 1994 Drop wavelength limits from calls to XCPLOT and XCRSLTS # Aug 12 1994 Change name of header-writing subrooutine # Aug 17 1994 Move ZPAD to labelled common # Aug 18 1994 Fix VCOMBINE call # Dec 8 1994 Move mwcs initialization to GETSPEC # Dec 19 1994 Read filter parameters here instead of in main program # Jan 31 1995 Change lengths of file and directory names to sz_pathname # Feb 15 1995 Print multi-pass results only if DEBUG is on # Feb 27 1995 Deal correctly with null spectrum and template directories # Mar 10 1995 Add option to read zero-padding flag from template # Mar 13 1995 Use CLOSE_IMAGE to make sure allocated space is freed # Mar 13 1995 Allow template-directed line chopping # Mar 13 1995 Rename zero_pad parameter for ease of relearning # Mar 15 1995 Allow template-directed template continuum removal # Mar 29 1995 Read rejection limits from xcsao.par, not contpars # Apr 4 1995 Move stack pointer setting to start of subroutine # Apr 6 1995 Change options for initial velocity setting # Apr 11 1995 Chop lines by default if template or spectrum header not set # May 8 1995 Don't add / to template directory if it ends in $ # May 15 1995 Change all sz_fname to sz_line, which is 100 chars longer # May 16 1995 Add template-driven option to divide, not subtract, continuum # May 18 1995 Zero best correlation before setting it # Jun 21 1995 Fix bug initializing velocity # Jun 26 1995 Add overlap header parameter to override wavelength limits # Jul 13 1995 Add debugging argument to VCOMBINE # Jul 19 1995 Call VWHEAD as well as XCWHEAD # Aug 11 1995 Set peak fraction from template header or parameter file # Sep 20 1995 Add CORRELATE parameter to bypass correlation for qplot # Oct 3 1995 If not CORRELATE, don't save results in header unless QPLOT # Oct 4 1995 Improve error message if no templates # Oct 13 1995 Always call SFREE before returning # Oct 13 1995 Print limits in angstroms; fix bug setting limits from .par # Feb 16 1996 Check filter limits to keep them within spectrum # Feb 22 1996 Set quality flag to 0 if correlating # Aug 7 1996 Use smw.h # Jan 31 1997 Add code to eliminate bad night sky lines # Feb 5 1997 If initial velocity is INDEF, start at zero # Mar 14 1997 Drop SPECSH from GETSPEC call; pass it in common # Mar 14 1997 Label spectrum dispersion axis from spectrum header # Apr 7 1997 Save ALL cross-correlations # Apr 14 1997 Add option to plot template spectra # Apr 25 1997 Add initial velocity option izguess for Z instead of velocity # Apr 30 1997 Pass NCOR to XCFILE so proper number of points is plotted # Aug 27 1997 Add band arguments for spectrum and template files # Dec 4 1997 Add spectrum and template aperture number to XCFILE call # Dec 22 1997 Make x-axis label work even if units aren't set in header # Feb 9 1998 Make maximum number of correlation points 16384, not 8192 # Feb 11 1998 Allocate correlation and velocity per template arrays # Feb 11 1998 Add code to count number of templates being correlated # May 15 1998 Allow minvel and maxvel to be INDEFD # Jun 12 1998 Use pixel limits of spectrum dispersion from header # Dec 29 1998 Fix bug when using full template pathnames in list # Mar 11 1999 Add write argument to getspec() # Mar 17 1999 Multiply renormalized fluxes by 1000