PRO plot_light_motion_curve, lmc_name, epoch, period, status, PLOT_EXP_MAGS = plot_exp_mags, PLOT_UNCAL = plot_uncal, PARALLAX = parallax ; Description: This module plots the light-motion curve "lmc_name" from the Light-Motion Curve Catalogue (LMCC - see ; Bramich et al. 2008). The separate light and motion curves are plotted in the postscript files ; ".Light.Curve.ps" and ".Motion.Curve.ps", respectively, which are created in the ; directory in which this module is run. ; For the light curve in each waveband, the module plots the time axis in units of days from the epoch ; "epoch", unless the parameter "period" is specified as a positive number, in which case the module plots ; the light curve in each waveband phased with the period "period" with zero phase at the epoch "epoch". ; Photometric data are plotted as black points with error bars. The mean magnitude of the light curve in ; each waveband is also plotted as a horizontal continuous line. ; For the motion curve along each coordinate axis (right ascension and declination), the module plots ; the time axis in units of days from the epoch "epoch". Astrometric data are plotted as black points with ; error bars. The astrometric data are fit with an astrometric model with constant proper motion. The fit ; for this model is plotted as a continuous line in the motion curve plots. ; This module provides the option to plot the light curve in each waveband using the exponential ; magnitudes instead of the PSF magnitudes, and it also provides the option to over-plot the uncalibrated ; magnitudes as red points with error bars. The module also provides the possibility to fit the motion ; curve with a second astrometric model that includes a constant proper motion and the term for parallax ; due to the Earth's orbital motion. When requested, the fit for this second model is over-plotted as a ; dashed line in the motion curve plots. ; ; Input Parameters: ; ; lmc_name - STRING - The name of the light-motion curve to be plotted, which should be of the form ; "lc_HH_MM_SS.SSS_pDD_MM_SS.SSS_t". The part of the name string with the format "HH_MM_SS.SSS" ; represents the right ascension (J2000) of the corresponding object in hours (HH, "00" to "23"), ; minutes (MM, "00" to "59"), and seconds (SS.SSS, "00.000" to "59.999"). The part of the name ; string with the format "pDD_MM_SS.SSS" represents the declination (J2000) of the corresponding ; object in degrees (DD, "00" to "90"), arcminutes (MM, "00" to "59"), and arcseconds (SS.SSS, ; "00.000" to "59.999") with the appropriate sign (p, "+" or "-"). ; epoch - FLOAT/DOUBLE - The heliocentric Julian date (days) to be used as zero on the time (or phase) axis in the ; light-motion curve plots. ; period - FLOAT/DOUBLE - The period (days) to be used to phase the light curve in each waveband. If this parameter ; is zero or negative, then the module will not phase the light curve. ; ; Output Parameters: ; ; status - INTEGER - If the module successfully plotted the light-motion curve, then "status" is returned with a ; value of "1", otherwise it is returned with a value of "0". ; ; Keywords: ; ; If the keyword PLOT_EXP_MAGS is set (as "/PLOT_EXP_MAGS"), then the module will plot the light curve in each ; waveband using the exponential magnitudes instead of the PSF magnitudes. ; ; If the keyword PLOT_UNCAL is set (as "/PLOT_UNCAL"), then the module will plot the light curve in each waveband ; using both the calibrated and uncalibrated magnitudes. The uncalibrated magnitudes are plotted with red points ; with error bars. This keyword applies to both the PSF and exponential magnitudes. ; ; If the keyword PARALLAX is set (as "/PARALLAX"), then the module will fit the motion curve with a second ; astrometric model that includes a constant proper motion and the term for parallax due to the Earth's orbital ; motion. The fit for this second model is over-plotted as a dashed line in the motion curve plots. ; ; Author: Dan Bramich (dan.bramich@hotmail.co.uk) ; ; Dependencies: (i) IDL Astronomy Users Library (http://idlastro.gsfc.nasa.gov/) ; (ii) DanIDL v2.0 (www.danidl.co.uk) ; ; History: ; ; 21/05/2013 - Module created (dmb). ;Perform basic checks on the input parameters print, '' print, 'Performing basic checks on the input parameters...' status = 0 if (test_str_scalar(lmc_name) NE 1) then begin print, 'ERROR - The parameter "lmc_name" is not a string...' return endif if (strlen(lmc_name) NE 31L) then begin print, 'ERROR - The parameter "lmc_name" is not a string of the correct length...' return endif if (test_fltdbl_scalar(epoch) NE 1) then begin print, 'ERROR - The parameter "epoch" is not a number of the correct type...' return endif if (test_fltdbl_scalar(period) NE 1) then begin print, 'ERROR - The parameter "period" is not a number of the correct type...' return endif epoch_use = double(epoch) period_use = double(period) ;Determine which file of the LMCC to read in file_str = 'LMCC.RA.' + strmid(lmc_name, 3, 2) + '.' + strmid(lmc_name, 6, 2) + '.fits' ;Read in the relevant LMCC file print, '' print, 'Reading in the LMCC file: ' + file_str if (file_test(file_str, /READ, /REGULAR) NE 1) then begin print, 'ERROR - The LMCC file "' + file_str + '" does not exist or it does not have the correct permissions...' return endif data = mrdfits(file_str, 1) ;Extract the data for the required light-motion curve print, '' print, 'Extracting the data for the light-motion curve: ' + lmc_name subs = where(data.light_curve_name EQ lmc_name, ndata) if (ndata EQ 0L) then begin print, 'ERROR - The light-motion curve "' + lmc_name + '" does not exist in the LMCC...' return endif if ((ndata mod 5L) NE 0L) then begin print, 'ERROR - The light-motion curve "' + lmc_name + '" is missing entries...' return endif data = data[subs] nepochs = ndata/5L filter = data.filter phot_hjd = data.hjd - epoch_use if (period_use GT 0.0D) then begin phot_phase = phot_hjd/period_use phot_phase = phot_phase - floor(phot_phase) endif print, 'Found ' + scmp(nepochs) + ' observational epochs...' ;Extract the PSF magnitudes or the exponential magnitudes, depending on which should be plotted print, 'Plotting the light curve for each waveband...' if keyword_set(plot_exp_mags) then begin mag_cal = double(data.exp_luptitude) mag_cal_err = double(data.exp_luptitude_error) endif else begin mag_cal = double(data.psf_luptitude) mag_cal_err = double(data.psf_luptitude_error) endelse ;Filter the calibrated photometric data for quality bit_patt1 = '20'X + '80'X + '40000'X + '80000'X + '400000'X bit_patt2 = '100'X + '800'X + '1000'X + '4000'X + '8000'X qual_tag = ((bit_patt1 AND data.flag1) EQ 0L) AND (('10000000'X AND data.flag1) NE 0L) AND ((bit_patt2 AND data.flag2) EQ 0L) good_cal_phot_tag = (data.photometric_calibration_tag EQ 1) AND (qual_tag EQ 1B) phot_cal_subs_u = where((filter EQ 0) AND (good_cal_phot_tag EQ 1B), nphot_cal_u) phot_cal_subs_g = where((filter EQ 1) AND (good_cal_phot_tag EQ 1B), nphot_cal_g) phot_cal_subs_r = where((filter EQ 2) AND (good_cal_phot_tag EQ 1B), nphot_cal_r) phot_cal_subs_i = where((filter EQ 3) AND (good_cal_phot_tag EQ 1B), nphot_cal_i) phot_cal_subs_z = where((filter EQ 4) AND (good_cal_phot_tag EQ 1B), nphot_cal_z) if ((nphot_cal_u + nphot_cal_g + nphot_cal_r + nphot_cal_i + nphot_cal_z) EQ 0L) then begin print, 'ERROR - The light-motion curve "' + lmc_name + '" does not contain any good data...' return endif ;Calculate the uncalibrated magnitudes by undoing the photometric calibration b_soft = double([1.4e-10, 0.9e-10, 1.2e-10, 1.8e-10, 7.4e-10]) mag_uncal = mag_cal mag_uncal_err = mag_cal_err subs = where(good_cal_phot_tag EQ 1B, nsubs) for i = 0L,(nsubs-1L) do begin csub = subs[i] curr_b_soft = b_soft[filter[csub]] result1 = lup_to_flux(mag_cal[csub], mag_cal_err[csub], curr_b_soft)*data[csub].photometric_zero_point result2 = flux_to_lup(result1[0], result1[1], curr_b_soft) mag_uncal[csub] = result2[0] mag_uncal_err[csub] = result2[1] endfor ;Filter the uncalibrated photometric data for quality phot_uncal_subs_u = where((filter EQ 0) AND (qual_tag EQ 1B), nphot_uncal_u) phot_uncal_subs_g = where((filter EQ 1) AND (qual_tag EQ 1B), nphot_uncal_g) phot_uncal_subs_r = where((filter EQ 2) AND (qual_tag EQ 1B), nphot_uncal_r) phot_uncal_subs_i = where((filter EQ 3) AND (qual_tag EQ 1B), nphot_uncal_i) phot_uncal_subs_z = where((filter EQ 4) AND (qual_tag EQ 1B), nphot_uncal_z) ;Determine the time-axis range hjd_lo = min(phot_hjd) hjd_hi = max(phot_hjd) hjd_range = (hjd_hi - hjd_lo) > 1.0D hjd_lo = hjd_lo - (0.1D*hjd_range) hjd_hi = hjd_hi + (0.1D*hjd_range) ;Determine the magnitude range mag_lo = dblarr(5) mag_hi = dblarr(5) if keyword_set(plot_uncal) then begin if (nphot_cal_u GT 0L) then begin mag_lo[0] = min([mag_cal[phot_cal_subs_u] - mag_cal_err[phot_cal_subs_u], mag_uncal[phot_uncal_subs_u] - mag_uncal_err[phot_uncal_subs_u]]) mag_hi[0] = max([mag_cal[phot_cal_subs_u] + mag_cal_err[phot_cal_subs_u], mag_uncal[phot_uncal_subs_u] + mag_uncal_err[phot_uncal_subs_u]]) endif else begin mag_lo[0] = min(mag_uncal[phot_uncal_subs_u] - mag_uncal_err[phot_uncal_subs_u]) mag_hi[0] = max(mag_uncal[phot_uncal_subs_u] + mag_uncal_err[phot_uncal_subs_u]) endelse if (nphot_cal_g GT 0L) then begin mag_lo[1] = min([mag_cal[phot_cal_subs_g] - mag_cal_err[phot_cal_subs_g], mag_uncal[phot_uncal_subs_g] - mag_uncal_err[phot_uncal_subs_g]]) mag_hi[1] = max([mag_cal[phot_cal_subs_g] + mag_cal_err[phot_cal_subs_g], mag_uncal[phot_uncal_subs_g] + mag_uncal_err[phot_uncal_subs_g]]) endif else begin mag_lo[1] = min(mag_uncal[phot_uncal_subs_g] - mag_uncal_err[phot_uncal_subs_g]) mag_hi[1] = max(mag_uncal[phot_uncal_subs_g] + mag_uncal_err[phot_uncal_subs_g]) endelse if (nphot_cal_r GT 0L) then begin mag_lo[2] = min([mag_cal[phot_cal_subs_r] - mag_cal_err[phot_cal_subs_r], mag_uncal[phot_uncal_subs_r] - mag_uncal_err[phot_uncal_subs_r]]) mag_hi[2] = max([mag_cal[phot_cal_subs_r] + mag_cal_err[phot_cal_subs_r], mag_uncal[phot_uncal_subs_r] + mag_uncal_err[phot_uncal_subs_r]]) endif else begin mag_lo[2] = min(mag_uncal[phot_uncal_subs_r] - mag_uncal_err[phot_uncal_subs_r]) mag_hi[2] = max(mag_uncal[phot_uncal_subs_r] + mag_uncal_err[phot_uncal_subs_r]) endelse if (nphot_cal_i GT 0L) then begin mag_lo[3] = min([mag_cal[phot_cal_subs_i] - mag_cal_err[phot_cal_subs_i], mag_uncal[phot_uncal_subs_i] - mag_uncal_err[phot_uncal_subs_i]]) mag_hi[3] = max([mag_cal[phot_cal_subs_i] + mag_cal_err[phot_cal_subs_i], mag_uncal[phot_uncal_subs_i] + mag_uncal_err[phot_uncal_subs_i]]) endif else begin mag_lo[3] = min(mag_uncal[phot_uncal_subs_i] - mag_uncal_err[phot_uncal_subs_i]) mag_hi[3] = max(mag_uncal[phot_uncal_subs_i] + mag_uncal_err[phot_uncal_subs_i]) endelse if (nphot_cal_z GT 0L) then begin mag_lo[4] = min([mag_cal[phot_cal_subs_z] - mag_cal_err[phot_cal_subs_z], mag_uncal[phot_uncal_subs_z] - mag_uncal_err[phot_uncal_subs_z]]) mag_hi[4] = max([mag_cal[phot_cal_subs_z] + mag_cal_err[phot_cal_subs_z], mag_uncal[phot_uncal_subs_z] + mag_uncal_err[phot_uncal_subs_z]]) endif else begin mag_lo[4] = min(mag_uncal[phot_uncal_subs_z] - mag_uncal_err[phot_uncal_subs_z]) mag_hi[4] = max(mag_uncal[phot_uncal_subs_z] + mag_uncal_err[phot_uncal_subs_z]) endelse endif else begin if (nphot_cal_u GT 0L) then begin mag_lo[0] = min(mag_cal[phot_cal_subs_u] - mag_cal_err[phot_cal_subs_u]) mag_hi[0] = max(mag_cal[phot_cal_subs_u] + mag_cal_err[phot_cal_subs_u]) endif if (nphot_cal_g GT 0L) then begin mag_lo[1] = min(mag_cal[phot_cal_subs_g] - mag_cal_err[phot_cal_subs_g]) mag_hi[1] = max(mag_cal[phot_cal_subs_g] + mag_cal_err[phot_cal_subs_g]) endif if (nphot_cal_r GT 0L) then begin mag_lo[2] = min(mag_cal[phot_cal_subs_r] - mag_cal_err[phot_cal_subs_r]) mag_hi[2] = max(mag_cal[phot_cal_subs_r] + mag_cal_err[phot_cal_subs_r]) endif if (nphot_cal_i GT 0L) then begin mag_lo[3] = min(mag_cal[phot_cal_subs_i] - mag_cal_err[phot_cal_subs_i]) mag_hi[3] = max(mag_cal[phot_cal_subs_i] + mag_cal_err[phot_cal_subs_i]) endif if (nphot_cal_z GT 0L) then begin mag_lo[4] = min(mag_cal[phot_cal_subs_z] - mag_cal_err[phot_cal_subs_z]) mag_hi[4] = max(mag_cal[phot_cal_subs_z] + mag_cal_err[phot_cal_subs_z]) endif endelse mag_range = max(mag_hi - mag_lo) mag_middle = 0.5D*(mag_lo + mag_hi) mag_lo = mag_middle - (0.55D*mag_range) mag_hi = mag_middle + (0.55D*mag_range) ;Set up the light curve plot set_plot, 'PS' device, XSIZE = 6.0, YSIZE = 10.0, /INCHES, /COLOR, BITS = 8, FILENAME = lmc_name + '.Light.Curve.ps' loadct, 12 dummy = dindgen(17)*(get_dbl_pi()/8.0D) usersym, 0.3D*cos(dummy), 0.3D*sin(dummy), /FILL !P.MULTI = [0, 1, 5] if (period_use GT 0.0D) then begin xrange = [-0.2, 1.2] xtitle = 'Phase' endif else begin xrange = [hjd_lo, hjd_hi] xtitle = 'HJD - ' + scmp(epoch_use, FORMAT = '(f25.1)') + ' (d)' endelse ;Plot the "u" waveband light curve plot, [0.0D], [0.0D], XRANGE = xrange, YRANGE = [mag_hi[0], mag_lo[0]], XSTYLE = 1, YSTYLE = 1, PSYM = 8, TITLE = lmc_name, $ XTITLE = xtitle, YTITLE = 'u', CHARSIZE = 1, /NODATA if (nphot_cal_u GT 0L) then begin if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_cal_subs_u], 0.0D], [mag_cal[phot_cal_subs_u], 0.0D], [mag_cal_err[phot_cal_subs_u], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_u] - 1.0D, 0.0D], [mag_cal[phot_cal_subs_u], 0.0D], [mag_cal_err[phot_cal_subs_u], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_u] + 1.0D, 0.0D], [mag_cal[phot_cal_subs_u], 0.0D], [mag_cal_err[phot_cal_subs_u], 0.0D], 8 endif else oploterr, [phot_hjd[phot_cal_subs_u], 0.0D], [mag_cal[phot_cal_subs_u], 0.0D], [mag_cal_err[phot_cal_subs_u], 0.0D], 8 result = weighted_arithmetic_mean(mag_cal[phot_cal_subs_u], mag_cal_err[phot_cal_subs_u], /ERRBAR) oplot, xrange, [result.mean, result.mean], LINESTYLE = 0 endif if keyword_set(plot_uncal) then begin !P.COLOR = 180 if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_uncal_subs_u], 0.0D], [mag_uncal[phot_uncal_subs_u], 0.0D], [mag_uncal_err[phot_uncal_subs_u], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_u] - 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_u], 0.0D], [mag_uncal_err[phot_uncal_subs_u], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_u] + 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_u], 0.0D], [mag_uncal_err[phot_uncal_subs_u], 0.0D], 8 endif else oploterr, [phot_hjd[phot_uncal_subs_u], 0.0D], [mag_uncal[phot_uncal_subs_u], 0.0D], [mag_uncal_err[phot_uncal_subs_u], 0.0D], 8 !P.COLOR = 0 endif ;Plot the "g" waveband light curve plot, [0.0D], [0.0D], XRANGE = xrange, YRANGE = [mag_hi[1], mag_lo[1]], XSTYLE = 1, YSTYLE = 1, PSYM = 8, TITLE = '', $ XTITLE = xtitle, YTITLE = 'g', CHARSIZE = 1, /NODATA if (nphot_cal_g GT 0L) then begin if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_cal_subs_g], 0.0D], [mag_cal[phot_cal_subs_g], 0.0D], [mag_cal_err[phot_cal_subs_g], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_g] - 1.0D, 0.0D], [mag_cal[phot_cal_subs_g], 0.0D], [mag_cal_err[phot_cal_subs_g], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_g] + 1.0D, 0.0D], [mag_cal[phot_cal_subs_g], 0.0D], [mag_cal_err[phot_cal_subs_g], 0.0D], 8 endif else oploterr, [phot_hjd[phot_cal_subs_g], 0.0D], [mag_cal[phot_cal_subs_g], 0.0D], [mag_cal_err[phot_cal_subs_g], 0.0D], 8 result = weighted_arithmetic_mean(mag_cal[phot_cal_subs_g], mag_cal_err[phot_cal_subs_g], /ERRBAR) oplot, xrange, [result.mean, result.mean], LINESTYLE = 0 endif if keyword_set(plot_uncal) then begin !P.COLOR = 180 if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_uncal_subs_g], 0.0D], [mag_uncal[phot_uncal_subs_g], 0.0D], [mag_uncal_err[phot_uncal_subs_g], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_g] - 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_g], 0.0D], [mag_uncal_err[phot_uncal_subs_g], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_g] + 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_g], 0.0D], [mag_uncal_err[phot_uncal_subs_g], 0.0D], 8 endif else oploterr, [phot_hjd[phot_uncal_subs_g], 0.0D], [mag_uncal[phot_uncal_subs_g], 0.0D], [mag_uncal_err[phot_uncal_subs_g], 0.0D], 8 !P.COLOR = 0 endif ;Plot the "r" waveband light curve plot, [0.0D], [0.0D], XRANGE = xrange, YRANGE = [mag_hi[2], mag_lo[2]], XSTYLE = 1, YSTYLE = 1, PSYM = 8, TITLE = '', $ XTITLE = xtitle, YTITLE = 'r', CHARSIZE = 1, /NODATA if (nphot_cal_r GT 0L) then begin if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_cal_subs_r], 0.0D], [mag_cal[phot_cal_subs_r], 0.0D], [mag_cal_err[phot_cal_subs_r], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_r] - 1.0D, 0.0D], [mag_cal[phot_cal_subs_r], 0.0D], [mag_cal_err[phot_cal_subs_r], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_r] + 1.0D, 0.0D], [mag_cal[phot_cal_subs_r], 0.0D], [mag_cal_err[phot_cal_subs_r], 0.0D], 8 endif else oploterr, [phot_hjd[phot_cal_subs_r], 0.0D], [mag_cal[phot_cal_subs_r], 0.0D], [mag_cal_err[phot_cal_subs_r], 0.0D], 8 result = weighted_arithmetic_mean(mag_cal[phot_cal_subs_r], mag_cal_err[phot_cal_subs_r], /ERRBAR) oplot, xrange, [result.mean, result.mean], LINESTYLE = 0 endif if keyword_set(plot_uncal) then begin !P.COLOR = 180 if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_uncal_subs_r], 0.0D], [mag_uncal[phot_uncal_subs_r], 0.0D], [mag_uncal_err[phot_uncal_subs_r], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_r] - 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_r], 0.0D], [mag_uncal_err[phot_uncal_subs_r], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_r] + 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_r], 0.0D], [mag_uncal_err[phot_uncal_subs_r], 0.0D], 8 endif else oploterr, [phot_hjd[phot_uncal_subs_r], 0.0D], [mag_uncal[phot_uncal_subs_r], 0.0D], [mag_uncal_err[phot_uncal_subs_r], 0.0D], 8 !P.COLOR = 0 endif ;Plot the "i" waveband light curve plot, [0.0D], [0.0D], XRANGE = xrange, YRANGE = [mag_hi[3], mag_lo[3]], XSTYLE = 1, YSTYLE = 1, PSYM = 8, TITLE = '', $ XTITLE = xtitle, YTITLE = 'i', CHARSIZE = 1, /NODATA if (nphot_cal_i GT 0L) then begin if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_cal_subs_i], 0.0D], [mag_cal[phot_cal_subs_i], 0.0D], [mag_cal_err[phot_cal_subs_i], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_i] - 1.0D, 0.0D], [mag_cal[phot_cal_subs_i], 0.0D], [mag_cal_err[phot_cal_subs_i], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_i] + 1.0D, 0.0D], [mag_cal[phot_cal_subs_i], 0.0D], [mag_cal_err[phot_cal_subs_i], 0.0D], 8 endif else oploterr, [phot_hjd[phot_cal_subs_i], 0.0D], [mag_cal[phot_cal_subs_i], 0.0D], [mag_cal_err[phot_cal_subs_i], 0.0D], 8 result = weighted_arithmetic_mean(mag_cal[phot_cal_subs_i], mag_cal_err[phot_cal_subs_i], /ERRBAR) oplot, xrange, [result.mean, result.mean], LINESTYLE = 0 endif if keyword_set(plot_uncal) then begin !P.COLOR = 180 if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_uncal_subs_i], 0.0D], [mag_uncal[phot_uncal_subs_i], 0.0D], [mag_uncal_err[phot_uncal_subs_i], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_i] - 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_i], 0.0D], [mag_uncal_err[phot_uncal_subs_i], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_i] + 1.0D, 0.0D], [mag_uncal[phot_uncal_subs_i], 0.0D], [mag_uncal_err[phot_uncal_subs_i], 0.0D], 8 endif else oploterr, [phot_hjd[phot_uncal_subs_i], 0.0D], [mag_uncal[phot_uncal_subs_i], 0.0D], [mag_uncal_err[phot_uncal_subs_i], 0.0D], 8 !P.COLOR = 0 endif ;Plot the "z" waveband light curve plot, [0.0D], [0.0D], XRANGE = xrange, YRANGE = [mag_hi[4], mag_lo[4]], XSTYLE = 1, YSTYLE = 1, PSYM = 8, TITLE = '', $ XTITLE = xtitle, YTITLE = 'z', CHARSIZE = 1, /NODATA if (nphot_cal_z GT 0L) then begin if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_cal_subs_z], 0.0D], [mag_cal[phot_cal_subs_z], 0.0D], [mag_cal_err[phot_cal_subs_z], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_z] - 1.0D, 0.0D], [mag_cal[phot_cal_subs_z], 0.0D], [mag_cal_err[phot_cal_subs_z], 0.0D], 8 oploterr, [phot_phase[phot_cal_subs_z] + 1.0D, 0.0D], [mag_cal[phot_cal_subs_z], 0.0D], [mag_cal_err[phot_cal_subs_z], 0.0D], 8 endif else oploterr, [phot_hjd[phot_cal_subs_z], 0.0D], [mag_cal[phot_cal_subs_z], 0.0D], [mag_cal_err[phot_cal_subs_z], 0.0D], 8 result = weighted_arithmetic_mean(mag_cal[phot_cal_subs_z], mag_cal_err[phot_cal_subs_z], /ERRBAR) oplot, xrange, [result.mean, result.mean], LINESTYLE = 0 endif if keyword_set(plot_uncal) then begin !P.COLOR = 180 if (period_use GT 0.0D) then begin oploterr, [phot_phase[phot_uncal_subs_z], 0.0D], [mag_uncal[phot_uncal_subs_z], 0.0D], [mag_uncal_err[phot_uncal_subs_z], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_z], 0.0D] - 1.0D, [mag_uncal[phot_uncal_subs_z], 0.0D], [mag_uncal_err[phot_uncal_subs_z], 0.0D], 8 oploterr, [phot_phase[phot_uncal_subs_z], 0.0D] + 1.0D, [mag_uncal[phot_uncal_subs_z], 0.0D], [mag_uncal_err[phot_uncal_subs_z], 0.0D], 8 endif else oploterr, [phot_hjd[phot_uncal_subs_z], 0.0D], [mag_uncal[phot_uncal_subs_z], 0.0D], [mag_uncal_err[phot_uncal_subs_z], 0.0D], 8 !P.COLOR = 0 endif ;Close the light curve plot file device, /CLOSE_FILE ;If there is only one observational epoch, then finish without plotting the motion curve if (nepochs EQ 1L) then begin print, 'WARNING - There is only one observational epoch and therefore the motion curve has not been plotted...' print, '' print, 'Finished!' status = 1 return endif ;Extract the astrometric measurements print, 'Plotting the motion curve for each coordinate axis...' subs = 5L*lindgen(nepochs) ast_hjd = phot_hjd[subs] for i = 1L,4L do begin subs = temporary(subs) + 1L ast_hjd = temporary(ast_hjd) + phot_hjd[subs] endfor ast_hjd = 0.2D*temporary(ast_hjd) ra = data[subs].ra dec = data[subs].dec ;Estimate the uncertainties on the astrometric measurements using the brightest uncalibrated magnitude with a good ;photometry flag over the five wavebands at each epoch brightest_mag_vec = replicate(10000.0D, nepochs) for i = 0L,(nepochs-1L) do begin csub = 5L*i for j = 0L,4L do begin if (qual_tag[csub + j] NE 1B) then continue brightest_mag_vec[i] = brightest_mag_vec[i] < mag_uncal[csub + j] endfor endfor ast_err_cut_hjd = 2453500.0D - epoch_use ra_dec_err = arcsec2deg(((0.0320D + (double(1.44e-14)*exp(1.34D*brightest_mag_vec)))*(ast_hjd LT ast_err_cut_hjd)) + $ ((0.0354D + (double(2.36e-12)*exp(1.09D*brightest_mag_vec)))*(ast_hjd GE ast_err_cut_hjd)), stat, /NO_PAR_CHECK) ra_dec_inv_var = 1.0D/(ra_dec_err^2) ra_dec_inv_var_sum = total(ra_dec_inv_var, /DOUBLE) ;Fit an astrometric model with constant proper motion to the astrometric data ra_dec_t0 = total(ast_hjd*ra_dec_inv_var, /DOUBLE)/ra_dec_inv_var_sum ra_dec_hjd = ast_hjd - ra_dec_t0 ra_dec_hjd_inv_var_sum = total((ra_dec_hjd*ra_dec_hjd)*ra_dec_inv_var, /DOUBLE) ra_mean = total(ra*ra_dec_inv_var, /DOUBLE)/ra_dec_inv_var_sum dec_mean = total(dec*ra_dec_inv_var, /DOUBLE)/ra_dec_inv_var_sum ra_pm = total(ra*ra_dec_hjd*ra_dec_inv_var, /DOUBLE)/ra_dec_hjd_inv_var_sum dec_pm = total(dec*ra_dec_hjd*ra_dec_inv_var, /DOUBLE)/ra_dec_hjd_inv_var_sum ;Correct the proper motion in right ascension for the declination ra_pm = ra_pm*cos(deg2rad(dec_mean, stat, /NO_PAR_CHECK)) ;If the keyword PARALLAX is set pl_model_tag = 0 if keyword_set(parallax) then begin ;If there are at least three observational epochs if (nepochs GT 2L) then begin ;Fit the astrometric data with an astrometric model that includes a constant proper motion and the term for ;parallax due to the Earth's orbital motion. The HJD 2451545.0 corresponds to the start of the year 2000. common data, d_ast_hjd, d_ra_dec_hjd, d_ra, d_dec, d_ra_dec_err d_ast_hjd = ast_hjd + (epoch_use - 2451545.0D) d_ra_dec_hjd = ra_dec_hjd d_ra = ra d_dec = dec d_ra_dec_err = ra_dec_err param = amoeba(1.0e-5, SCALE = [5.0D/3600.0D, 5.0D/3600.0D, 5.0D/(365.25D*3600.0D), 5.0D/(365.25D*3600.0D), (1.0D/3600.0D)], $ P0 = [ra_mean, dec_mean, ra_pm, dec_pm, 0.0D], FUNCTION_NAME = 'chisqr_par') pl_ra_mean = param[0] pl_dec_mean = param[1] pl_ra_pm = param[2] pl_dec_pm = param[3] pl = param[4] pl_model_tag = 1 ;If there are only two observational epochs, then issue a warning that it was not possible to fit the second ;astrometric model that includes a constant proper motion and the term for parallax due to the Earth's orbital ;motion endif else print, 'WARNING - Cannot fit an astrometric model with constant proper motion and parallax due to a lack of observational epochs...' endif ;Determine the right ascension and declination range ra = deg2arcsec(ra - ra_mean, stat, /NO_PAR_CHECK) ra_dec_err = deg2arcsec(ra_dec_err, stat, /NO_PAR_CHECK) ra_lo = min(ra - ra_dec_err) ra_hi = max(ra + ra_dec_err) ra_range = (ra_hi - ra_lo) > 0.01D ra_lo = ra_lo - (0.1D*ra_range) ra_hi = ra_hi + (0.1D*ra_range) dec = deg2arcsec(dec - dec_mean, stat, /NO_PAR_CHECK) dec_lo = min(dec - ra_dec_err) dec_hi = max(dec + ra_dec_err) dec_range = (dec_hi - dec_lo) > 0.01D dec_lo = dec_lo - (0.1D*dec_range) dec_hi = dec_hi + (0.1D*dec_range) ra_dec_lo = min([ra_lo, dec_lo]) ra_dec_hi = max([ra_hi, dec_hi]) ;Generate the fitted astrometric models for plotting grid_hjd = ((dindgen(1001)/1000.0D)*(hjd_hi - hjd_lo)) + hjd_lo ra_mean = deg2rad(ra_mean, stat, /NO_PAR_CHECK) dec_mean = deg2rad(dec_mean, stat, /NO_PAR_CHECK) ra_pm = deg2arcsec(ra_pm, stat, /NO_PAR_CHECK) dec_pm = deg2arcsec(dec_pm, stat, /NO_PAR_CHECK) ra_model = (ra_pm/cos(dec_mean))*(grid_hjd - ra_dec_t0) dec_model = dec_pm*(grid_hjd - ra_dec_t0) if (pl_model_tag EQ 1) then begin L = rev2rad(((grid_hjd + (epoch_use - 2451545.0D))/365.25D) + 0.75D, stat, /NO_PAR_CHECK) pl_ra_mean = deg2rad(pl_ra_mean, stat, /NO_PAR_CHECK) pl_dec_mean = deg2rad(pl_dec_mean, stat, /NO_PAR_CHECK) pl_ra_pm = deg2arcsec(pl_ra_pm, stat, /NO_PAR_CHECK) pl_dec_pm = deg2arcsec(pl_dec_pm, stat, /NO_PAR_CHECK) epsilon = deg2rad(23.439281D, stat, /NO_PAR_CHECK) pl = deg2arcsec(pl, stat, /NO_PAR_CHECK) ra_par_fac = ((cos(epsilon)*cos(pl_ra_mean)*sin(L)) - (sin(pl_ra_mean)*cos(L)))/cos(pl_dec_mean) pl_ra_model = rad2arcsec(pl_ra_mean - ra_mean, stat, /NO_PAR_CHECK) + ((pl_ra_pm/cos(pl_dec_mean))*(grid_hjd - ra_dec_t0)) + (pl*ra_par_fac) dec_par_fac = (sin(epsilon)*cos(pl_dec_mean)*sin(L)) - (sin(pl_dec_mean)*cos(pl_ra_mean)*cos(L)) - (cos(epsilon)*sin(pl_dec_mean)*sin(pl_ra_mean)*sin(L)) pl_dec_model = rad2arcsec(pl_dec_mean - dec_mean, stat, /NO_PAR_CHECK) + (pl_dec_pm*(grid_hjd - ra_dec_t0)) + (pl*dec_par_fac) endif ;Set up the motion curve plot set_plot, 'PS' device, XSIZE = 6.0, YSIZE = 4.0, /INCHES, /COLOR, BITS = 8, FILENAME = lmc_name + '.Motion.Curve.ps' !P.MULTI = [0, 1, 2] pm_str = scmp(365.25D*sqrt((ra_pm^2) + (dec_pm^2)), FORMAT='(f10.6)') ;Plot the right ascension motion curve plot, [0.0D], [0.0D], XRANGE = [hjd_lo, hjd_hi], YRANGE = [ra_dec_lo, ra_dec_hi], XSTYLE = 1, YSTYLE = 1, PSYM = 8, $ TITLE = lmc_name, XTITLE = 'HJD - ' + scmp(epoch_use, FORMAT = '(f25.1)') + ' (d)', YTITLE = 'RA - Mean RA (arcsec)', $ CHARSIZE = 0.6, /NODATA oploterr, ast_hjd, ra, ra_dec_err, 8 oplot, grid_hjd, ra_model, LINESTYLE = 0 if (pl_model_tag EQ 1) then oplot, grid_hjd, pl_ra_model, LINESTYLE = 1 ;Plot the declination motion curve plot, [0.0D], [0.0D], XRANGE = [hjd_lo, hjd_hi], YRANGE = [ra_dec_lo, ra_dec_hi], XSTYLE = 1, YSTYLE = 1, PSYM = 8, $ TITLE = 'Proper Motion (arcsec/yr): ' + pm_str, XTITLE = 'HJD - ' + scmp(epoch_use, FORMAT = '(f25.1)') + ' (d)', $ YTITLE = 'Dec - Mean Dec (arcsec)', CHARSIZE = 0.6, /NODATA oploterr, ast_hjd, dec, ra_dec_err, 8 oplot, grid_hjd, dec_model, LINESTYLE = 0 if (pl_model_tag EQ 1) then oplot, grid_hjd, pl_dec_model, LINESTYLE = 1 ;Close the motion curve plot file device, /CLOSE_FILE ;Set "status" to "1" print, '' print, 'Finished!' status = 1 END