diff -Nru psrchive-2012.12+20170930/Base/Classes/Profile.C psrchive-2012.12+20171106/Base/Classes/Profile.C --- psrchive-2012.12+20170930/Base/Classes/Profile.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Classes/Profile.C 2017-11-06 13:33:09.000000000 +0000 @@ -377,7 +377,7 @@ float log_threshold = log(threshold)/log(base); - if (!finite(log_threshold)) + if (!isfinite(log_threshold)) throw Error (InvalidParam, "Pulsar::Profile::logarithm", "logarithm of threshold=%lf is not finite", threshold); @@ -390,7 +390,7 @@ amps[ibin] = log(amps[ibin])/log(base); else amps[ibin] = log_threshold; - if (!finite(amps[ibin])) + if (!isfinite(amps[ibin])) throw Error (InvalidParam, "Pulsar::Profile::logarithm", "logarithm of amps[%u]=%lf is not finite", ibin, amps[ibin]); diff -Nru psrchive-2012.12+20170930/Base/Classes/Profile_rotate.C psrchive-2012.12+20171106/Base/Classes/Profile_rotate.C --- psrchive-2012.12+20170930/Base/Classes/Profile_rotate.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Classes/Profile_rotate.C 2017-11-06 13:33:09.000000000 +0000 @@ -33,7 +33,7 @@ */ void Pulsar::Profile::rotate_phase (double phase) { - if (!finite(phase)) + if (!isfinite(phase)) throw Error (InvalidParam, "Pulsar::Profile::rotate_phase", "non-finite phase = %lf\n", phase); diff -Nru psrchive-2012.12+20170930/Base/Classes/Pulsar/Integration.h psrchive-2012.12+20171106/Base/Classes/Pulsar/Integration.h --- psrchive-2012.12+20170930/Base/Classes/Pulsar/Integration.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Classes/Pulsar/Integration.h 2017-11-06 13:33:09.000000000 +0000 @@ -509,6 +509,22 @@ for (unsigned ichan=0; ichanget_Profile(ipol, ichan)->*(method)) (arg); } + + template + void foreach (Integration* integration, const Integration* operand, + BinaryProfileMethod method) + { + const unsigned npol = integration->get_npol(); + const unsigned nchan = integration->get_nchan(); + + for (unsigned ipol=0; ipolget_Profile(ipol, ichan); + const Profile* from = operand->get_Profile(ipol, ichan); + (into->*(method)) (from); + } + } } #endif diff -Nru psrchive-2012.12+20170930/Base/Formats/EPN/EPNArchive.C psrchive-2012.12+20171106/Base/Formats/EPN/EPNArchive.C --- psrchive-2012.12+20170930/Base/Formats/EPN/EPNArchive.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Formats/EPN/EPNArchive.C 2017-11-06 13:33:09.000000000 +0000 @@ -368,8 +368,8 @@ "error calling crwpen"); if (verbose == 3) { - cerr << "Pulsar::EPNArchive::read_record rwepn called" << endl; - epn_dump (&line1, &line2, &line3, &line4, &line5); + cerr << "Pulsar::EPNArchive::read_record crwepn called" << endl; + // epn_dump (&line1, &line2, &line3, &line4, &line5); } current_record = record; diff -Nru psrchive-2012.12+20170930/Base/Formats/EPN/epnio.c psrchive-2012.12+20171106/Base/Formats/EPN/epnio.c --- psrchive-2012.12+20170930/Base/Formats/EPN/epnio.c 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Formats/EPN/epnio.c 2017-11-06 13:33:09.000000000 +0000 @@ -207,20 +207,6 @@ return 0; } -void epn_dump (const epn_header_line1* line1, - const epn_header_line2* line2, - const epn_header_line3* line3, - const epn_header_line4* line4, - const epn_header_line5* line5) -{ - fprintf (stderr, "line 5: npol=%d\n", line5->npol); - fprintf (stderr, "line 5: nfreq=%d\n", line5->nfreq); - fprintf (stderr, "line 5: nbin=%d\n", line5->nbin); - fprintf (stderr, "line 5: tbin=%lf\n", line5->tbin); - fprintf (stderr, "line 5: nint=%d\n", line5->nint); - fprintf (stderr, "line 5: ncal=%d\n", line5->ncal); -} - int cnepnrec (const char* filename) { /* for dealing with filename string */ diff -Nru psrchive-2012.12+20170930/Base/Formats/Timer/fcomp.C psrchive-2012.12+20171106/Base/Formats/Timer/fcomp.C --- psrchive-2012.12+20170930/Base/Formats/Timer/fcomp.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Formats/Timer/fcomp.C 2017-11-06 13:33:09.000000000 +0000 @@ -55,7 +55,7 @@ N_FromLittleEndian (nvals, packed_buf); } - if (scale==0 || !finite(scale)) + if (scale==0 || !isfinite(scale)) { cerr << "fcompread: invalid scale=" << scale << endl; return -1; diff -Nru psrchive-2012.12+20170930/Base/Formats/Timer/TimerArchive_load.C psrchive-2012.12+20171106/Base/Formats/Timer/TimerArchive_load.C --- psrchive-2012.12+20170930/Base/Formats/Timer/TimerArchive_load.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Formats/Timer/TimerArchive_load.C 2017-11-06 13:33:09.000000000 +0000 @@ -239,7 +239,7 @@ if (hdr.corrected == 1) hdr.corrected = NINT_CORRECTED; - if (hdr.corrected & NINT_CORRECTED == 0) + if ((hdr.corrected & NINT_CORRECTED) == 0) { /* nint() correction jss */ double difftime,tcorr; @@ -272,7 +272,7 @@ if (hdr.corrected == 1) hdr.corrected = S2ROS_CORRECTED; - if (hdr.corrected & S2ROS_CORRECTED == 0) + if ((hdr.corrected & S2ROS_CORRECTED) == 0) { /* get the year data was processed */ char comment[64]; diff -Nru psrchive-2012.12+20170930/Base/Formats/Timer/TimerProfile_unload.C psrchive-2012.12+20171106/Base/Formats/Timer/TimerProfile_unload.C --- psrchive-2012.12+20170930/Base/Formats/Timer/TimerProfile_unload.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Formats/Timer/TimerProfile_unload.C 2017-11-06 13:33:09.000000000 +0000 @@ -34,7 +34,7 @@ const float* amps = profile->get_amps(); for (int ibin=0; ibin < nbin; ibin++) - if (!finite(amps[ibin])) + if (!isfinite(amps[ibin])) throw Error (InvalidParam, "Pulsar::TimerProfile_unload", "amps[%d]=%f is not finite", ibin, amps[ibin]); diff -Nru psrchive-2012.12+20170930/Base/Formats/WAPP/WAPPArchive.C psrchive-2012.12+20171106/Base/Formats/WAPP/WAPPArchive.C --- psrchive-2012.12+20170930/Base/Formats/WAPP/WAPPArchive.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Base/Formats/WAPP/WAPPArchive.C 2017-11-06 13:33:09.000000000 +0000 @@ -203,11 +203,12 @@ throw Error (FailedSys, "Pulsar::WAPPArchive::load_header", "fseek(EOF)"); } + + if (ftell(f) < 0) + throw Error (FailedSys, "Pulsar::WAPPArchive::load_header", "ftell"); + wapp_file_size = ftell(f); fclose(f); - if (wapp_file_size<0) - throw Error (FailedSys, "Pulsar::WAPPArchive::load_header", - "ftell()"); // Check that this is a folding-mode file. if (strncmp(hdr->obs_type, "PULSAR_FOLDING", 15)!=0) diff -Nru psrchive-2012.12+20170930/debian/changelog psrchive-2012.12+20171106/debian/changelog --- psrchive-2012.12+20170930/debian/changelog 2017-10-23 14:43:39.000000000 +0000 +++ psrchive-2012.12+20171106/debian/changelog 2017-11-08 08:25:47.000000000 +0000 @@ -1,3 +1,9 @@ +psrchive (2012.12+20171106-1kern1) xenial; urgency=medium + + * Imported Upstream version 2012.12+20171106 + + -- KERN packaging Wed, 08 Nov 2017 10:25:27 +0200 + psrchive (2012.12+20170930-1kern1) xenial; urgency=medium * rebuild diff -Nru psrchive-2012.12+20170930/More/Applications/pat.C psrchive-2012.12+20171106/More/Applications/pat.C --- psrchive-2012.12+20170930/More/Applications/pat.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Applications/pat.C 2017-11-06 13:33:09.000000000 +0000 @@ -84,7 +84,7 @@ #if HAVE_PGPLOT void plotDifferences(Pulsar::Archive* arch, Pulsar::Archive* stdarch, vector& toas, const double min_phase, const double max_phase, - const bool output_plot_difference); + const bool output_profile_residuals); void set_phase_zoom(vector >& plots, const double min, const double max); @@ -191,6 +191,7 @@ " -r Print reference phase and dt \n" " -R Print only the phase shift and error in turns \n" " -u Print as pat-like format smjd + dt \n" + " -b Output profile residuals \n" "\n" "Plotting options (if compiled with pgplot):\n" " -K Specify plot device\n" @@ -227,10 +228,10 @@ bool skip_bad = false; bool phase_info = false; bool tempo2_output = false; + bool output_profile_residuals = false; #if HAVE_PGPLOT bool plot_difference = false; - bool output_plot_difference = false; bool centre_template_peak = false; #endif @@ -302,11 +303,9 @@ break; } -#if HAVE_PGPLOT case 'b': - output_plot_difference = true; + output_profile_residuals = true; break; -#endif case 'c': choose_maximum_harmonic = true; @@ -625,6 +624,9 @@ arrival->set_observation (arch); + if (output_profile_residuals) + arrival->set_residual( arch->clone() ); + if (verbose) cerr << "pat: calling ArrivalTime::get_toas" << endl; @@ -639,10 +641,17 @@ arch->remove_baseline(); rotate_archive(arch, toas); plotDifferences(arch, stdarch, toas, min_phase, max_phase, - output_plot_difference); + output_profile_residuals); } #endif + if (output_profile_residuals) + { + Archive* residual = arrival->get_residual (); + string filename = arch->get_filename() + ".resid"; + residual->unload (filename); + } + if (phase_only) { for (unsigned i = 0; i < toas.size(); i++) @@ -786,14 +795,14 @@ * @param stdarch Template archive. * @param min_phase min x-value when zooming * @param max_phase max x-value when zooming - * @param output_plot_difference whether or not the difference profile + * @param output_profile_residuals whether or not the difference profile * should be written out to .diff */ #if HAVE_PGPLOT void plotDifferences(Pulsar::Archive* arch, Pulsar::Archive* stdarch, vector& toas, const double min_phase, const double max_phase, - const bool output_plot_difference) + const bool output_profile_residuals) { // remove baseline for all templates (except caldelay) @@ -814,7 +823,7 @@ plotter->configure("x:range=" + get_xrange(min_phase, max_phase)); ofstream f; - if (output_plot_difference) { + if (output_profile_residuals) { string output_filename = replace_extension(arch->get_filename(), plot_difference_extension); @@ -877,7 +886,7 @@ profile_plot(plotter, profile_archive, profile, arch, centre_frequency); - if (output_plot_difference) { + if (output_profile_residuals) { // Output (to .diff): // float* bins = profile_diff->get_Profile(0,0,0)->get_amps(); diff -Nru psrchive-2012.12+20170930/More/Applications/paz.C psrchive-2012.12+20171106/More/Applications/paz.C --- psrchive-2012.12+20170930/More/Applications/paz.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Applications/paz.C 2017-11-06 13:33:09.000000000 +0000 @@ -376,7 +376,8 @@ while (fgets (buffer, 4096, fptr)) { char* key = strtok (buffer, whitespace); - while (key) + // parse tokens until a comment is started + while (key && key[0] != '#') { chans_to_zero.push_back (fromstring(key)); key = strtok (NULL, whitespace); diff -Nru psrchive-2012.12+20170930/More/Applications/psrzap.C psrchive-2012.12+20171106/More/Applications/psrzap.C --- psrchive-2012.12+20170930/More/Applications/psrzap.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Applications/psrzap.C 2017-11-06 13:33:09.000000000 +0000 @@ -193,9 +193,10 @@ << endl << "Mouse:" << endl << " Left-click selects the start of a range" << endl - << " then left-click again to zoom, or right-click to zap." << endl - << " Right-click zaps current cursor location." << endl - << " Middle-click (or 'd') updates the diagnostic plots." << endl + << " then left-click again to zoom, or right-click (or 'X') to zap." + << endl + << " Right-click (or 'X') zaps current cursor location." << endl + << " Middle-click (or 'D') updates the diagnostic plots." << endl << "Keyboard:" << endl << " " << CMD_HELP << " Show this help" << endl << " " << CMD_FREQMODE << " Use frequency select mode" << endl diff -Nru psrchive-2012.12+20170930/More/Applications/rmfit.C psrchive-2012.12+20171106/More/Applications/rmfit.C --- psrchive-2012.12+20170930/More/Applications/rmfit.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Applications/rmfit.C 2017-11-06 13:33:09.000000000 +0000 @@ -1091,9 +1091,9 @@ { double dm = data->get_dispersion_measure(); float dm_auto_maxrm = auto_maxrm_dm * dm; - cerr << "rmfit: maximum RM = (" << auto_maxrm_dm << " * DM=" << dm - << ") = " << dm_auto_maxrm << endl; - auto_maxrm = std::min (auto_maxrm, dm_auto_maxrm); + cerr << "rmfit: maximum RM = (" << auto_maxrm_dm << " * DM=" << dm + << ") = " << dm_auto_maxrm << endl; + auto_maxrm = std::min (auto_maxrm, dm_auto_maxrm); } bool override = false; @@ -1110,8 +1110,8 @@ { cerr << "rmfit: default auto maximum RM under limit = " << auto_minrm << endl; - RM_max = auto_minrm; - override = true; + RM_max = auto_minrm; + override = true; } if (auto_step_rad > 0) @@ -1133,23 +1133,23 @@ override = true; } - if (override) - { - rmsteps = unsigned( 2 * RM_max / step_size ); - cerr << "rmfit: override maximum RM=" << RM_max - << " steps=" << rmsteps << endl; - } - - if (auto_minsteps > 0 && rmsteps < auto_minsteps) - { - rmsteps = auto_minsteps; - cerr << "rmfit: override steps=" << auto_minsteps - << " (step size=" << 2*RM_max/rmsteps << ")" << endl; - } + if (override) + { + rmsteps = unsigned( 2 * RM_max / step_size ); + cerr << "rmfit: override maximum RM=" << RM_max + << " steps=" << rmsteps << endl; + } - minrm = -RM_max; - maxrm = RM_max; + if (auto_minsteps > 0 && rmsteps < auto_minsteps) + { + rmsteps = auto_minsteps; + cerr << "rmfit: override steps=" << auto_minsteps + << " (step size=" << 2*RM_max/rmsteps << ")" << endl; } + + minrm = -RM_max; + maxrm = RM_max; + } vector fluxes (rmsteps); vector rms (rmsteps); @@ -1159,37 +1159,51 @@ Reference::To backup = data->clone(); + double max_snr = 0.0; + double max_L = 0.0; + for (unsigned step=0; step < rmsteps; step++) - { - double rm = minrm + step * rmstepsize; - - /* - WvS, 12 July 2006: Note that Archive::defaraday can handle the - case that the data have already been corrected for a certain - rotation measure. - - Wvs, 26 September 2007: Then again, perhaps round-off error - can build up over many iterations. - */ - - data->set_rotation_measure( rm ); - data->defaraday (); + { + double rm = minrm + step * rmstepsize; - Reference::To useful = data->clone(); - useful->fscrunch(); - useful->remove_baseline(); + /* + WvS, 12 July 2006: Note that Archive::defaraday can handle the + case that the data have already been corrected for a certain + rotation measure. + + Wvs, 26 September 2007: Then again, perhaps round-off error + can build up over many iterations. + */ + + data->set_rotation_measure( rm ); + data->defaraday (); - poln_stats.set_profile( useful->get_Integration(0)->new_PolnProfile(0) ); - Estimate rval = poln_stats.get_total_linear (); + Reference::To useful = data->clone(); + useful->fscrunch(); + useful->remove_baseline(); + + Reference::To profile; + profile = useful->get_Integration(0)->new_PolnProfile(0); + + poln_stats.set_profile( profile ); + Estimate rval = poln_stats.get_total_linear (); + + fluxes[step] = rval.get_value(); + err[step] = rval.get_error(); + rms[step] = rm; - fluxes[step] = rval.get_value(); - err[step] = rval.get_error(); - rms[step] = rm; + if (fluxes[step] > max_L) + { + max_L = fluxes[step]; - data = backup->clone(); + Pulsar::Profile linear; + profile->get_linear (&linear); + max_snr = linear.snr(); } - - + + data = backup->clone(); + } + ofstream os ("rm_spectrum.txt"); for (unsigned step=0; step < rmsteps; step++) os << rms[step] << " " << fluxes[step] << " " << err[step] << endl; @@ -1325,8 +1339,14 @@ cerr << "Width="<< gm.get_width() <<" Height="<< gm.get_height() << endl; + double FWHM = 2.0 * sqrt ( 2.0 * log (2.0) ) * gm.get_width(); + cerr << "rmfit: FWHM/SNR uncertainty=" << 0.5 * FWHM / max_snr << endl; + bestrm = gm.get_centre(); + Estimate estimatedRM = gm.get_Estimate (0); + cerr << "rmfit: estimated RM = " << estimatedRM << endl; + // Compute the error based on how far you have to move the // centre of the Gaussian before the two-sigma diff -Nru psrchive-2012.12+20170930/More/General/Statistics.C psrchive-2012.12+20171106/More/General/Statistics.C --- psrchive-2012.12+20170930/More/General/Statistics.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/General/Statistics.C 2017-11-06 13:33:09.000000000 +0000 @@ -185,32 +185,48 @@ } //! Get the off-pulse baseline -Pulsar::PhaseWeight* Pulsar::Statistics::get_baseline () +Pulsar::PhaseWeight* Pulsar::Statistics::get_baseline () try { setup_stats (); return stats->get_baseline(); } + catch (Error& error) + { + throw error += "Pulsar::Statistics::get_baseline"; + } //! Get the on-pulse phase bins -Pulsar::PhaseWeight* Pulsar::Statistics::get_onpulse () +Pulsar::PhaseWeight* Pulsar::Statistics::get_onpulse () try { setup_stats (); return stats->get_onpulse(); } + catch (Error& error) + { + throw error += "Pulsar::Statistics::get_onpulse"; + } //! Get all phase bins -Pulsar::PhaseWeight* Pulsar::Statistics::get_all () +Pulsar::PhaseWeight* Pulsar::Statistics::get_all () try { setup_stats (); return stats->get_all(); } + catch (Error& error) + { + throw error += "Pulsar::Statistics::get_all"; + } //! Get the statistics calculator -Pulsar::ProfileStats* Pulsar::Statistics::get_stats () +Pulsar::ProfileStats* Pulsar::Statistics::get_stats () try { setup_stats (); return stats; } + catch (Error& error) + { + throw error += "Pulsar::Statistics::get_stats"; + } //! Get the weighted frequency of the Pulsar::Archive double Pulsar::Statistics::get_weighted_frequency () const @@ -244,7 +260,7 @@ return dispersion_smear (dm, freq, chan_bw); } -void Pulsar::Statistics::setup_stats () +void Pulsar::Statistics::setup_stats () try { if (!stats) { @@ -252,6 +268,8 @@ stats_setup = false; } + // avoid recursion - part 2 + // (Plugin::setup might call a function that calls setup_stats) if (stats_setup) return; @@ -263,7 +281,8 @@ if (Profile::verbose) cerr << "Pulsar::Statistics::setup_stats profile=" << profile.ptr() << endl; - // avoid recursion (Plugin::setup might call a function that calls setup_stats) + // avoid recursion - part 1 + // (Plugin::setup might call a function that calls setup_stats) stats_setup = true; for (unsigned i=0; i 2) - cerr << "HybridCalibrator response=" << response << endl; + cerr << "HybridCalibrator pre-calibrator response=" << response << endl; } // pass the reference Stokes parameters through the instrument @@ -389,15 +389,22 @@ solver->set_output (output_stokes); solver->solve (correction); - // produce the supplemented transformation, + MEAL::Complex2* result = correction; - MEAL::ProductRule* result; - result = new MEAL::ProductRule; + if (precalibrator) + { + /* the result is the product of the correction + and the precalibrator response */ - result->add_model( correction ); - if (precalibrator) - result->add_model( new MEAL::Value(response) ); + MEAL::ProductRule* product; + product = new MEAL::ProductRule; + product->add_model( correction ); + product->add_model( new MEAL::Value(response) ); + + result = product; + } + transformation[ichan] = result; } diff -Nru psrchive-2012.12+20170930/More/Polarimetry/Pulsar/ReceptionModelSolveCeres.h psrchive-2012.12+20171106/More/Polarimetry/Pulsar/ReceptionModelSolveCeres.h --- psrchive-2012.12+20170930/More/Polarimetry/Pulsar/ReceptionModelSolveCeres.h 1970-01-01 00:00:00.000000000 +0000 +++ psrchive-2012.12+20171106/More/Polarimetry/Pulsar/ReceptionModelSolveCeres.h 2017-11-06 13:33:09.000000000 +0000 @@ -0,0 +1,38 @@ +//-*-C++-*- +/*************************************************************************** + * + * Copyright (C) 2017 by Willem van Straten + * Licensed under the Academic Free License version 2.1 + * + ***************************************************************************/ + +// psrchive/More/Polarimetry/Pulsar/ReceptionModelSolveCeres.h + +#ifndef __ReceptionModel_SolveCeres_H +#define __ReceptionModel_SolveCeres_H + +#include "Pulsar/ReceptionModelSolver.h" + +namespace Calibration +{ + //! Solve the measurement equation using Google's Ceres Solver + class SolveCeres : public ReceptionModel::Solver + { + + //! Return the name of this solver + std::string get_name () const; + + //! Return a new, copy-constructed clone + SolveCeres* clone () const; + + protected: + + //! Solve the measurement equation using Ceres Solver + void fit (); + + }; + +} + +#endif + diff -Nru psrchive-2012.12+20170930/More/Polarimetry/ReceptionModelSolveCeres.C psrchive-2012.12+20171106/More/Polarimetry/ReceptionModelSolveCeres.C --- psrchive-2012.12+20170930/More/Polarimetry/ReceptionModelSolveCeres.C 1970-01-01 00:00:00.000000000 +0000 +++ psrchive-2012.12+20171106/More/Polarimetry/ReceptionModelSolveCeres.C 2017-11-06 13:33:09.000000000 +0000 @@ -0,0 +1,39 @@ +/*************************************************************************** + * + * Copyright (C) 2017 by Willem van Straten + * Licensed under the Academic Free License version 2.1 + * + ***************************************************************************/ + +#include "Pulsar/ReceptionModelSolveCeres.h" +#include "Pulsar/CoherencyMeasurementSet.h" + +#include "ceres/ceres.h" + +#include +#include + +using namespace std; +using Calibration::CoherencyMeasurementSet; + +std::string Calibration::SolveCeres::get_name () const +{ + return "Ceres"; +} + +Calibration::SolveCeres* Calibration::SolveCeres::clone () const +{ + return new SolveCeres; +} + +void Calibration::SolveCeres::fit () +{ + if (Calibration::ReceptionModel::verbose) + cerr << "Calibration::SolveCeres::Fit" + " nfit=" << get_nparam_infit () << + " ndat=" << get_ndat_constraint () << endl; + + + ceres::Problem problem; +} + diff -Nru psrchive-2012.12+20170930/More/Timing/ArrivalTime.C psrchive-2012.12+20171106/More/Timing/ArrivalTime.C --- psrchive-2012.12+20170930/More/Timing/ArrivalTime.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Timing/ArrivalTime.C 2017-11-06 13:33:09.000000000 +0000 @@ -10,7 +10,7 @@ #include "Pulsar/PolnProfileShiftEstimator.h" #include "Pulsar/Archive.h" -#include "Pulsar/Integration.h" +#include "Pulsar/IntegrationExpert.h" #include "Pulsar/Profile.h" #include "Pulsar/IntegrationOrder.h" #include "Pulsar/Statistics.h" @@ -163,7 +163,8 @@ const Profile* profile = subint->get_Profile (0, ichan); - if (skip_bad && (profile->get_weight() == 0) || (shift && standard->get_Profile (0,0,ichan)->get_weight() == 0)) + if ((skip_bad && (profile->get_weight() == 0)) + || (shift && standard->get_Profile (0,0,ichan)->get_weight() == 0)) continue; try @@ -183,6 +184,15 @@ arrival_time.set_arrival(arrival_time.get_arrival() + be->get_delay()); toas.push_back( arrival_time ); + + if (residual) + { + Integration* rsubint = residual->get_Integration (isub); + rsubint->expert()->rotate_phase( shift.get_value() ); + + const Integration* std = standard->get_Integration (0); + foreach (rsubint, std, &Profile::diff); + } } catch (Error& error) { diff -Nru psrchive-2012.12+20170930/More/Timing/MatrixTemplateMatching.C psrchive-2012.12+20171106/More/Timing/MatrixTemplateMatching.C --- psrchive-2012.12+20170930/More/Timing/MatrixTemplateMatching.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Timing/MatrixTemplateMatching.C 2017-11-06 13:33:09.000000000 +0000 @@ -1,6 +1,6 @@ /*************************************************************************** * - * Copyright (C) 2009 by Willem van Straten + * Copyright (C) 2009 - 2017 by Willem van Straten * Licensed under the Academic Free License version 2.1 * ***************************************************************************/ @@ -11,6 +11,7 @@ #include "Pulsar/Archive.h" #include "Pulsar/Integration.h" +#include "Pulsar/PolnProfile.h" using namespace std; @@ -64,7 +65,7 @@ std::vector& toas) { engine->add_pulsar (observation, isub); - + const Integration* integration = observation->get_Integration (isub); unsigned nchan = integration->get_nchan(); @@ -73,12 +74,28 @@ if (!engine->get_transformation_valid(ichan)) continue; - Estimate shift = engine->get_mtm(ichan) -> get_phase(); + const PolnProfileFit* mtm = engine->get_mtm(ichan); + + Estimate shift = mtm -> get_phase(); Tempo::toa TOA = get_toa (shift, integration, ichan); - TOA.set_reduced_chisq( engine->get_mtm(ichan)->get_reduced_chisq () ); + TOA.set_reduced_chisq( mtm->get_reduced_chisq () ); toas.push_back( TOA ); + + if (residual) + { + Jones xform = mtm->get_transformation()->evaluate(); + + Integration* rsubint = residual->get_Integration (isub); + Reference::To rprof = rsubint->new_PolnProfile (ichan); + rprof->transform( inv(xform) ); + rprof->rotate_phase( shift.get_value() ); + + const Integration* std = standard->get_Integration (0); + Reference::To stdprof = std->new_PolnProfile (ichan); + rprof->diff (stdprof); + } } #if 0 diff -Nru psrchive-2012.12+20170930/More/Timing/Pulsar/ArrivalTime.h psrchive-2012.12+20171106/More/Timing/Pulsar/ArrivalTime.h --- psrchive-2012.12+20170930/More/Timing/Pulsar/ArrivalTime.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Timing/Pulsar/ArrivalTime.h 2017-11-06 13:33:09.000000000 +0000 @@ -77,10 +77,23 @@ //! Get auxilliary information std::string get_value (const std::string& key, const Tempo::toa&); + //! Set the archive that will store residual profiles + void set_residual (Archive* res) { residual = res; } + + //! Get the archive of residual profiles + Archive* get_residual () { return residual; } + protected: + //! The observation to be fit to the standard Reference::To observation; + + //! The standard to which observations are fit Reference::To standard; + + //! The residual pulse profiles (transformed observation minus standard) + Reference::To residual; + Reference::To shift_estimator; //! default TOA output format diff -Nru psrchive-2012.12+20170930/More/Timing/ShiftEstimator.C psrchive-2012.12+20171106/More/Timing/ShiftEstimator.C --- psrchive-2012.12+20170930/More/Timing/ShiftEstimator.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/More/Timing/ShiftEstimator.C 2017-11-06 13:33:09.000000000 +0000 @@ -44,5 +44,8 @@ /*! Most estimators work with total intensity */ void Pulsar::ShiftEstimator::preprocess (Archive* archive) { - archive->pscrunch (); + if (archive->get_npol() == 4) + archive->convert_state( Signal::Stokes ); + else + archive->pscrunch (); } diff -Nru psrchive-2012.12+20170930/Util/fft/QuaternionFT.C psrchive-2012.12+20171106/Util/fft/QuaternionFT.C --- psrchive-2012.12+20170930/Util/fft/QuaternionFT.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/fft/QuaternionFT.C 2017-11-06 13:33:09.000000000 +0000 @@ -231,10 +231,10 @@ out += 4; } - delete H_a; - delete H_b; - delete h_a; - delete h_b; + delete [] H_a; + delete [] H_b; + delete [] h_a; + delete [] h_b; } Quaternion exp (const Quaternion& axis, float phi) diff -Nru psrchive-2012.12+20170930/Util/fitsutil/psrfitsio.C psrchive-2012.12+20171106/Util/fitsutil/psrfitsio.C --- psrchive-2012.12+20170930/Util/fitsutil/psrfitsio.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/fitsutil/psrfitsio.C 2017-11-06 13:33:09.000000000 +0000 @@ -13,6 +13,8 @@ using namespace std; +bool psrfits_verbose = false; + static void update_tdim (fitsfile* ffptr, int column, unsigned ndim, ...) { vector dims (ndim); diff -Nru psrchive-2012.12+20170930/Util/fitsutil/psrfitsio.h psrchive-2012.12+20171106/Util/fitsutil/psrfitsio.h --- psrchive-2012.12+20170930/Util/fitsutil/psrfitsio.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/fitsutil/psrfitsio.h 2017-11-06 13:33:09.000000000 +0000 @@ -19,6 +19,8 @@ #include #include +extern bool psrfits_verbose; + //! Empty template class requires specialization template struct FITS_traits { }; @@ -197,12 +199,20 @@ int colnum = 0; int status = 0; + if (psrfits_verbose) + std::cerr << "psrfits_write_col calling fits_get_colnum" + " name='" << name << "'" << std::endl; + fits_get_colnum (fptr, CASEINSEN, const_cast(name), &colnum, &status); if (status) throw FITSError (status, "psrfits_write_col(vector)", "fits_get_colnum (name=%s)", name); + if (psrfits_verbose) + std::cerr << "psrfits_write_col calling fits_modify_vector_len" + " colnum=" << colnum << " size=" << data.size() << std::endl; + fits_modify_vector_len (fptr, colnum, data.size(), &status); if (status) @@ -210,9 +220,15 @@ "fits_modify_vector_len (name=%s col=%d size=%u)", name, colnum, data.size()); + if (psrfits_verbose) + std::cerr << "psrfits_write_col calling psrfits_update_tdim" << std::endl; + if (dims.size() > 1) psrfits_update_tdim (fptr, colnum, dims); + if (psrfits_verbose) + std::cerr << "psrfits_write_col calling fits_write_col" << std::endl; + fits_write_col (fptr, FITS_traits::datatype(), colnum, row, 1, data.size(), diff -Nru psrchive-2012.12+20170930/Util/genutil/angleconv.c psrchive-2012.12+20171106/Util/genutil/angleconv.c --- psrchive-2012.12+20170930/Util/genutil/angleconv.c 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/angleconv.c 2017-11-06 13:33:09.000000000 +0000 @@ -137,38 +137,6 @@ *sec = radt < 0.0 ? 0.0 : (radt * 60) ; } -void makebname( double * ra, double * raerr, double * dec, double * decerr, char * name) -{ - double sec, tmp = *ra/15; - int deg, min; - char isign; - - name[0]='B'; - rad2dms(&tmp, &isign, °, &min, &sec); - sprintf(name+1,"%02d%02d",deg,min); - rad2dms(dec, &isign, °, &min, &sec); - sprintf(name+5,"%c%02d",isign,deg); - name[8]='\0'; -} - - -void makejname( double * ra, double * raerr, double * dec, double * decerr, char * name) -{ - double sec, tmp = *ra/15; - int deg, min; - char isign; - - name[0]='J'; - rad2dms(&tmp, &isign, °, &min, &sec); - sprintf(name+1,"%02d%02d",deg,min); - rad2dms(dec, &isign, °, &min, &sec); - sprintf(name+5,"%c%02d%02d",isign,deg,min); - if(*decerr > 0.0 && (*decerr / M_PI * 180.0 * 60) > 1.0 ) /* truncate */ - fprintf(stderr,"Recommend renaming to %s\n",name); /* name if necessary */ - /* name[8]='\0'; */ - name[10]='\0'; -} - char * crad2dms( double pos, int ra, char * name ) { double sec, tmp; diff -Nru psrchive-2012.12+20170930/Util/genutil/dwt_undec.c psrchive-2012.12+20171106/Util/genutil/dwt_undec.c --- psrchive-2012.12+20170930/Util/genutil/dwt_undec.c 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/dwt_undec.c 2017-11-06 13:33:09.000000000 +0000 @@ -31,7 +31,7 @@ void udwt_forward_step(double *in, double *outl, double *outh, size_t n, const double h[], const double g[], size_t nc, int s) { - int i,j,idx; + size_t i,j,idx; /* Zero output array */ for (i=0; i=0; i--) { + for (int i=nl-1; i>=0; i--) { udwt_inverse_step(curil, curih, curo, n, h, g, nc, i); swp = curil; curil = curo; diff -Nru psrchive-2012.12+20170930/Util/genutil/f77util.c psrchive-2012.12+20171106/Util/genutil/f77util.c --- psrchive-2012.12+20170930/Util/genutil/f77util.c 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/f77util.c 2017-11-06 13:33:09.000000000 +0000 @@ -6,7 +6,9 @@ ***************************************************************************/ #include "f77util.h" + #include +#include void f2cstr (const char* f_str, char* c_str, unsigned length) { diff -Nru psrchive-2012.12+20170930/Util/genutil/ThreadMemory.C psrchive-2012.12+20171106/Util/genutil/ThreadMemory.C --- psrchive-2012.12+20170930/Util/genutil/ThreadMemory.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/ThreadMemory.C 2017-11-06 13:33:09.000000000 +0000 @@ -21,25 +21,39 @@ ThreadMemory::ThreadMemory () { +#ifdef HAVE_PTHREAD pthread_key_create (&key, &destructor); +#else + memory = 0; +#endif } ThreadMemory::~ThreadMemory () { +#ifdef HAVE_PTHREAD pthread_key_delete (key); +#else + destructor (memory); +#endif } void * ThreadMemory::get () { - +#ifdef HAVE_PTHREAD void * memory= pthread_getspecific (key); if (!memory) throw Error (InvalidState, "ThreadMemory::get", "key ptr not set for this thread"); +#endif return memory; } -void ThreadMemory::set (void * memory) +void ThreadMemory::set (void* _memory) { - pthread_setspecific (key, memory); +#ifdef HAVE_PTHREAD + pthread_setspecific (key, _memory); +#else + memory = _memory; +#endif } + diff -Nru psrchive-2012.12+20170930/Util/genutil/ThreadMemory.h psrchive-2012.12+20171106/Util/genutil/ThreadMemory.h --- psrchive-2012.12+20170930/Util/genutil/ThreadMemory.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/ThreadMemory.h 2017-11-06 13:33:09.000000000 +0000 @@ -40,6 +40,8 @@ #ifdef HAVE_PTHREAD pthread_key_t key; +#else + void* memory; #endif }; diff -Nru psrchive-2012.12+20170930/Util/genutil/ThreadStream.C psrchive-2012.12+20171106/Util/genutil/ThreadStream.C --- psrchive-2012.12+20170930/Util/genutil/ThreadStream.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/ThreadStream.C 2017-11-06 13:33:09.000000000 +0000 @@ -20,26 +20,40 @@ ThreadStream::ThreadStream () { +#ifdef HAVE_PTHREAD pthread_key_create (&key, &destructor); +#else + stream = 0; +#endif } ThreadStream::~ThreadStream () { +#ifdef HAVE_PTHREAD pthread_key_delete (key); +#else + if (stream) delete stream; +#endif } ostream& ThreadStream::get () { - ostream* value = reinterpret_cast( pthread_getspecific (key) ); +#ifdef HAVE_PTHREAD + ostream* stream = reinterpret_cast( pthread_getspecific (key) ); - if (!value) + if (!stream) throw Error (InvalidState, "ThreadStream::get", "std::ostream not set for this thread"); - return *value; +#endif + return *stream; } -void ThreadStream::set (std::ostream* stream) +void ThreadStream::set (std::ostream* _stream) { - pthread_setspecific (key, stream); +#ifdef HAVE_PTHREAD + pthread_setspecific (key, _stream); +#else + stream = _stream; +#endif } diff -Nru psrchive-2012.12+20170930/Util/genutil/ThreadStream.h psrchive-2012.12+20171106/Util/genutil/ThreadStream.h --- psrchive-2012.12+20170930/Util/genutil/ThreadStream.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/genutil/ThreadStream.h 2017-11-06 13:33:09.000000000 +0000 @@ -48,6 +48,8 @@ #ifdef HAVE_PTHREAD pthread_key_t key; +#else + std::ostream* stream; #endif }; diff -Nru psrchive-2012.12+20170930/Util/tempo/inverse_phase.h psrchive-2012.12+20171106/Util/tempo/inverse_phase.h --- psrchive-2012.12+20170930/Util/tempo/inverse_phase.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/tempo/inverse_phase.h 2017-11-06 13:33:09.000000000 +0000 @@ -62,7 +62,7 @@ dt = (predictor.phase(guess) - p) / predictor.frequency(guess); - if (!finite(dt.in_seconds())) + if (!isfinite(dt.in_seconds())) throw Error (InvalidState, "inverse_phase", "dt not finite; freq=%lf", predictor.frequency(guess)); diff -Nru psrchive-2012.12+20170930/Util/tempo/Phase.C psrchive-2012.12+20171106/Util/tempo/Phase.C --- psrchive-2012.12+20170930/Util/tempo/Phase.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/tempo/Phase.C 2017-11-06 13:33:09.000000000 +0000 @@ -83,7 +83,7 @@ string s = tostring (turns); string f = tostring (fturns, precision, std::ios::fixed); - if (!finite(fturns)) + if (!isfinite(fturns)) s += ".NaN"; else { diff -Nru psrchive-2012.12+20170930/Util/tempo2/T2Predictor.C psrchive-2012.12+20171106/Util/tempo2/T2Predictor.C --- psrchive-2012.12+20170930/Util/tempo2/T2Predictor.C 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/tempo2/T2Predictor.C 2017-11-06 13:33:09.000000000 +0000 @@ -254,7 +254,7 @@ "epoch %s not spanned by ChebyModelSet", t.printdays(20).c_str()); - if (!finite(p)) { + if (!isfinite(p)) { Error error (InvalidState, "Tempo2::Predictor::phase", "T2Predictor_GetPhase result = "); error << p; diff -Nru psrchive-2012.12+20170930/Util/third/gnu/dummy.c psrchive-2012.12+20171106/Util/third/gnu/dummy.c --- psrchive-2012.12+20170930/Util/third/gnu/dummy.c 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/third/gnu/dummy.c 1970-01-01 00:00:00.000000000 +0000 @@ -1,7 +0,0 @@ -/* This file is only a place holder to placate ar on OS X */ - -static int do_nothing () -{ - -} - diff -Nru psrchive-2012.12+20170930/Util/third/gnu/Makefile.am psrchive-2012.12+20171106/Util/third/gnu/Makefile.am --- psrchive-2012.12+20170930/Util/third/gnu/Makefile.am 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/third/gnu/Makefile.am 2017-11-06 13:33:09.000000000 +0000 @@ -1,8 +1,8 @@ noinst_LTLIBRARIES = libgnu.la -libgnu_la_SOURCES = dummy.c include_HEADERS = +libgnu_la_SOURCES = if !HAVE_GETOPT_LONG libgnu_la_SOURCES += getopt.c getopt1.c diff -Nru psrchive-2012.12+20170930/Util/third/Makefile.am psrchive-2012.12+20171106/Util/third/Makefile.am --- psrchive-2012.12+20170930/Util/third/Makefile.am 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/third/Makefile.am 2017-11-06 13:33:09.000000000 +0000 @@ -1,5 +1,5 @@ -SUBDIRS = parsifal star gnu +SUBDIRS = parsifal star noinst_LTLIBRARIES = libthird.la @@ -7,8 +7,12 @@ libthird_la_LIBADD = \ parsifal/libparsifal.la \ - star/libstar.la \ - gnu/libgnu.la + star/libstar.la + +if !HAVE_GETOPT_LONG +SUBDIRS += gnu +libthird_la_LIBADD += gnu/libgnu.la +endif include $(top_srcdir)/config/Makefile.local diff -Nru psrchive-2012.12+20170930/Util/units/tostring.h psrchive-2012.12+20171106/Util/units/tostring.h --- psrchive-2012.12+20170930/Util/units/tostring.h 2017-10-23 13:07:28.000000000 +0000 +++ psrchive-2012.12+20171106/Util/units/tostring.h 2017-11-06 13:33:09.000000000 +0000 @@ -23,6 +23,8 @@ { typedef std::ios_base::fmtflags fmtflags; +#define fmtzero fmtflags(0) + unsigned precision; bool precision_set; @@ -36,7 +38,12 @@ ToString () { reset_modifiers(); } - void reset_modifiers () { precision_set = setf_set = unsetf_set = false; } + void reset_modifiers () + { + precision_set = setf_set = unsetf_set = false; + precision = 0; + setf = unsetf = fmtzero; + } void set_precision (unsigned p) { precision = p; precision_set = true; } void set_setf (fmtflags f) { setf = f; setf_set = true; }