text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include "psrsalsa.h"
typedef struct {
  datafile_definition datafile;
  double f2_min, f2_max;
  double f3_min, f3_max;
  double P3IntegrateLow, P3IntegrateHigh;
  int SelectP3Integrate;
  double oversaturize;
}twodfs_def;
typedef struct {
  char *textside;
  double fl_min, fl_max;
  double dl;
  double LRegionLow, LRegionHigh;
  int SelectLRegion;
  int doFlip;
  int ExtraVerticalMaximaSkip;
  int noside;
  int overlaypp;
  int inside;
  int noylabels;
  int noxlabels;
  int notop;
  int autooversaturizel, autooversaturize2, autozoomP2;
  float labelscale, titlescale;
  double oversaturizel;
  double Imin, Imax;
  double scaleFig_x, scaleFig_y;
  int plot_xlabel, plot_ylabel, plot_ylabeltop;
  int nomain;
  int showwedge;
  int showwedge_max;
  int normaliseSide;
  int nointegrateNumbers;
  int intflip;
  int usephase;
  int lineColour;
  int normalise_spectra;
  int SelectP3Region;
}plotoptions_def;
int loadLRFS(datafile_definition *lrfs, int extprefix, int longsnap, int argc, char **argv, plotoptions_def *plotoptions, verbose_definition verbose);
int load2dfs(twodfs_def *twodfs_allinfo, datafile_definition AverageProfile, int file_number, int extprefix, int altname, int argc, char **argv, plotoptions_def plotoptions, verbose_definition verbose);
void Plot2dfs(twodfs_def twodfs_allinfo, twodfs_def twodfs2_allinfo, datafile_definition AverageProfile, int twodfsonly, int Number, char *title, plotoptions_def plotoptions, verbose_definition verbose);
void PlotLRFS(datafile_definition lrfs, datafile_definition AverageProfile, twodfs_def twodfs_allinfo, char *title, plotoptions_def plotoptions, verbose_definition verbose);
int loadHeaderPulseStack(datafile_definition *AverageProfile, datafile_definition lrfs, double *period, int type_of_plots, int argc, char **argv, verbose_definition verbose);
void autozoomP2(twodfs_def *twodfs_allinfo, datafile_definition AverageProfile, plotoptions_def plotoptions, verbose_definition verbose);
int main(int argc, char **argv)
{
  int i, j, xi, SelectP2Region, LoadTwo, NrSelectedOverSaturize;
  int file_number, ImaxSet, IminSet, type_of_plots, maxSubpulsePhaseSet, minSubpulsePhaseSet, do_phase_slope, extprefix;
  char PlotDevice[100], filename[1000], txt[1000];
  int title_index, Load2dfs, LoadLRFS, plotvariance, plotmodindex, ok_flag, altProf, lineStyle, ret, longsnap;
  double P3RegionLow, P3RegionHigh, maxSubpulsePhase, minSubpulsePhase;
  float I, x, profilescale;
  float maxvalue_mod, maxvalue_stddev, maxsigma_stddev, maxsigma_mod, ImaxValue, IminValue, phase_slope_g, phase_slope_o;
  FILE *fout_ascii;
  long k;
  datafile_definition subpulseTrackProfile, subpulseTrackProfileErr, subpulseAmpProfile;
  datafile_definition AverageProfile, VarianceProfile, ModProfile, VarianceProfileErr, ModProfileErr, lrfs;
  verbose_definition noverbose;
  psrsalsaApplication application;
  plotoptions_def plotoptions;
  twodfs_def twodfs_allinfo, twodfs2_allinfo;
  initApplication(&application, "pspecFig", "");
  type_of_plots = 0;
  file_number = 1;
  plotoptions.scaleFig_x = 1;
  plotoptions.scaleFig_y = 1;
  do_phase_slope = 0;
  plotoptions.autooversaturizel = 0;
  plotoptions.autooversaturize2 = 0;
  plotoptions.autozoomP2 = 0;
  extprefix = 0;
  longsnap = 0;
  if(argc < 2) {
    printf("Program to plot some of the pspec output\n\nUsage: pspecFig [options] stack_file (i.e. the name of the pulse stack that has been processed by pspec). By default the profile/modulation index/standard deviation profile/lrfs/2dfs are combined in a single plot (mode A). When -phaseplot is specified, a plot of the profile/subpulse amplitude and subpulse phase is produced (mode B).\n\n");
    printf("Where optional options are:\n\n");
    printf("General options:\n");
    printf("-headerlist       Show options of the -header option.\n");
    printf("-header           Change this header parameter.\n");
    printf("                  Example -header 'name J0123+4567'\n");
    printf("-v                Verbose mode (to get a better idea what is happening)\n");
    printf("-debug            Enable more output (where implemented)\n");
    printf("-nocounters       Don't show counters etc (useful when generating log files)\n");
    printf("-phaseplot        Enable mode B operation, see above.\n");
    printf("\nGeneral file/panel options:\n");
    printf("-altProf          Load this alternative .profile file (generated by pspec).\n");
    printf("                  Useful if profile/modindex/stddev curves are calculated from\n");
    printf("                  a different file than the spectra.\n");
    printf("-notop            Don't show top plot\n");
    printf("\nMode A specific file/panel options:\n");
    printf("-2                Load two 2dfs's corresponding to two pulse longitude ranges\n");
    printf("-2dfsnr           Override 2dfs file number of first 2dfs. Default is %d.\n", file_number);
    printf("-nolrfs           Do not load lrfs\n");
    printf("-no2dfs           Do not load 2dfs's\n");
    printf("-noside           Don't show side panels\n");
    printf("-normspectra      Normalise the spectra (peak value=1)\n");
    printf("-normside         Normalise the side panels (peak value=1)\n");
    printf("\nGeneral range options:\n");
    printf("-Imax             Set maximum value of the y-range of the top plot\n");
    printf("-Imin             Set minimum value of the y-range of the top plot\n");
    printf("-long  \"low high\" Set horizontal range shown lrfs/profile.\n");
    printf("-phase            Use pulse longitude in phase rather than degrees.\n");
    printf("-dlong            Shift lrfs/profile by this amount of degrees or phase.\n");
    printf("\nMode A specific range options:\n");
    printf("-p3 \"low high\"    Set vertical range shown in lrfs/2dfs in cpp.\n");
    printf("-p2 \"low high\"    Set horizontal range as shown in 2dfs (can use this option\n");
    printf("                  twice if -2 is used) in cpp.\n");
    printf("-int  \"low high\"  Select vertical integration range in 2dfs in cpp\n");
    printf("                  (affects the side panels, can be use twice if -2 is used).\n");
    printf("-modsigma         Set minimum significance for the modulation index (def. is 3).\n");
    printf("-stddevsigma      Same for the stddev values (default is 3).\n");
    printf("-modmax           Set the max allowed value for the modulation index\n");
    printf("                  (other values are ignored).\n");
    printf("-stddevmax        Same for the stddev values\n");
    printf("\nMode B specific range options:\n");
    printf("-spmax            Set maximum value of the y-range of the subpulse phase plot\n");
    printf("-spmin            Set minimum value of the y-range of the subpulse phase plot\n");
    printf("\nOther graphics options:\n");
    printf("-device  \"...\"    Plot device\n");
    printf("-inside           Place tick marks inside\n");
    printf("-labelscale       Set size of labels (default is 1)\n");
    printf("-titlescale       Set size of title (default is 1), but -labalscale affects size as well.\n");
    printf("-noylabels        Don't show ylabels\n");
    printf("-linestyle        Set the PGPLOT line style of the pulse profile\n");
    printf("-linecolor        Set the PGPLOT line color of the pulse profile\n");
    printf("-scalefig \"x y\"   Scale size of panel with factors x and y\n");
    printf("-title            Set the title\n");
    printf("-ytop             Show y-label of top plot\n");
    printf("\nMode A specific graphics options:\n");
    printf("-intflip          Add a flipped version of the vertical integration of the 2dfs\n");
    printf("-intnrs           Show nrs along side panels axis, instead of just a tick\n");
    printf("-scalel scale     The values in the lrfs are boosted by factor scale, resulting\n");
    printf("                  in clipping, which can highlight weaker features.\n");
    printf("-scale2 scale     As -scalel, but for 2dfs instead of lrfs.\n");
    printf("                  When this option is used twice, the 2nd time the\n");
    printf("                  option is used applies to the second 2dfs shown.\n");
    printf("-scalep           Scale profile by this factor\n");
    printf("-noflip           Do not flip 2DFS horizontally. If specified positive drift\n");
    printf("                  corresponds to power in the negative side of the diagram.\n");
    printf("-overlay          Overlay pulse profile over LRFS\n");
    printf("-noxlabels        Don't show xlabels on 2dfs bottom integration panel.\n");
    printf("-nomod            Do not plot a modulation index profile\n");
    printf("-nostddev         Do not plot a standard deviation profile\n");
    printf("-showwedge        Plot an annotated wedge to show color scale\n");
    printf("-showwedge_showmax  Plot a label indicating the maximum value in the data.\n");
    printf("                  With the -scalel or -scale2 option this is not necessarily\n");
    printf("                  the maximum of the colour scale.\n");
    printf("-textside         Set the text printed in top left corner of the graph\n");
    printf("-xlabel           Show x-label of LRFS and 2DFS in cpp\n");
    printf("-xlabel2          Show x-label of LRFS and 2DFS in P0/P2\n");
    printf("-xlabel3          Show x-label of LRFS and 2DFS in P/P2\n");
    printf("-ylabel           Show y-label of LRFS in cpp\n");
    printf("-ylabel2          Show y-label of LRFS in P0/P3\n");
    printf("-ylabel3          Show y-label of LRFS in P/P3\n");
    printf("-s                Skip the specified number of P3 bins when determining\n");
    printf("                  the range in left-hand side integrations panels. By default\n");
    printf("                  the first bin is skipped.\n");
    printf("\nMode B specific graphics options:\n");
    printf("-noflip           Do not change sign of subpulse phase. If specified positive\n");
    printf("                  drift corresponds to a decrease subpulse phase as function of\n");
    printf("                  pulse longitude.\n");
    printf("-xlabel           Show x-label in pulse longitude\n");
    printf("-ylabel           Show y-label in subpulse phase\n");
    printf("-phaseslope \"g o\" In the subpulse phase plot, add a line with gradient g (in deg\n");
    printf("                  per deg) and offset o in deg.\n");
    printf("\n");
    printf("Please use the appropriate citation when using results of this software in your publications:\n\n");
    printf("More information about the lrfs/2dfs/modulation index can be found in:\n");
    printf(" - Weltevrede et al. 2006, A&A, 445, 243\n");
    printf(" - Weltevrede et al. 2007, A&A, 469, 607\n");
    printf("More information about bootstrap/subpulse phase track & amplitude can be found in:\n");
    printf(" - Weltevrede et al. 2012, MNRAS, 424, 843\n\n");
    printCitationInfo();
    terminateApplication(&application);
    return 0;
  }
  plotoptions.plot_ylabeltop = 0;
  plotoptions.nomain = 0;
  plotoptions.SelectP3Region = 0;
  SelectP2Region = 0;
  plotoptions.SelectLRegion = 0;
  twodfs_allinfo.SelectP3Integrate = 0;
  twodfs2_allinfo.SelectP3Integrate = 0;
  LoadTwo = 0;
  strcpy(PlotDevice, "?");
  twodfs_allinfo.oversaturize = 1;
  twodfs2_allinfo.oversaturize = 1;
  plotoptions.oversaturizel = 1;
  NrSelectedOverSaturize = 0;
  maxvalue_mod = -1;
  maxsigma_mod = 3;
  maxvalue_stddev = -1;
  maxsigma_stddev = 3;
  title_index = -1;
  plotoptions.textside = NULL;
  plotoptions.doFlip = 1;
  plotoptions.ExtraVerticalMaximaSkip = 1;
  Load2dfs = 1;
  LoadLRFS = 1;
  plotvariance = 1;
  plotmodindex = 1;
  profilescale = 1;
  plotoptions.plot_xlabel = 0;
  plotoptions.plot_ylabel = 0;
  plotoptions.notop = 0;
  plotoptions.noside = 0;
  plotoptions.overlaypp = 0;
  plotoptions.inside = 0;
  plotoptions.labelscale = 1.0;
  plotoptions.titlescale = 1.0;
  plotoptions.noylabels = 0;
  plotoptions.noxlabels = 0;
  altProf = 0;
  plotoptions.dl = 0;
  ImaxSet = 0;
  IminSet = 0;
  lineStyle = 1;
  plotoptions.lineColour = 1;
  maxSubpulsePhaseSet = 0;
  minSubpulsePhaseSet = 0;
  maxSubpulsePhase = 0;
  minSubpulsePhase = 0;
  plotoptions.showwedge = 0;
  plotoptions.showwedge_max = 0;
  plotoptions.normalise_spectra = 0;
  plotoptions.normaliseSide = 0;
  plotoptions.nointegrateNumbers = 1;
  plotoptions.intflip = 0;
  plotoptions.usephase = 0;
  twodfs_allinfo.f2_min = 0.0;
  twodfs_allinfo.f2_max = 0.0;
  twodfs2_allinfo.f2_min = 0.0;
  twodfs2_allinfo.f2_max = 0.0;
  for(i = 1; i < argc; i++) {
    if(strcmp(argv[i], "-headerlist") == 0) {
      printHeaderCommandlineOptions(stdout);
      terminateApplication(&application);
      return 0;
    }
  }
  int argclast;
  argclast = argc-1;
  if(argv[argc-1][0] == '-')
    argclast += 1;
  for(i = 1; i < argclast; i++) {
    if(strcmp(argv[i], "-d") == 0 || strcmp(argv[i], "-D") == 0 || strcmp(argv[i], "-device") == 0 || strcmp(argv[i], "-dev") == 0) {
      strcpy(PlotDevice, argv[i+1]);
      i++;
    }else if(strcmp(argv[i], "-v") == 0) {
      application.verbose_state.verbose = 1;
    }else if(strcmp(argv[i], "-phaseplot") == 0) {
      type_of_plots = 1;
    }else if(strcmp(argv[i], "-debug") == 0) {
      application.verbose_state.debug = 1;
    }else if(strcmp(argv[i], "-nocounters") == 0) {
      application.verbose_state.nocounters = 1;
    }else if(strcmp(argv[i], "-xlabel") == 0
      ) {
      plotoptions.plot_xlabel = 1;
    }else if(strcmp(argv[i], "-xlabel2") == 0
      ) {
      plotoptions.plot_xlabel = 2;
    }else if(strcmp(argv[i], "-xlabel3") == 0
      ) {
      plotoptions.plot_xlabel = 3;
    }else if(strcmp(argv[i], "-ylabel") == 0
      ) {
      plotoptions.plot_ylabel = 1;
    }else if(strcmp(argv[i], "-ylabel2") == 0
      ) {
      plotoptions.plot_ylabel = 2;
    }else if(strcmp(argv[i], "-ylabel3") == 0
      ) {
      plotoptions.plot_ylabel = 3;
    }else if(strcmp(argv[i], "-nolrfs") == 0
      ) {
      LoadLRFS = 0;
    }else if(strcmp(argv[i], "-2") == 0) {
      LoadTwo = 1;
    }else if(strcmp(argv[i], "-no2dfs") == 0
      ) {
      Load2dfs = 0;
    }else if(strcmp(argv[i], "-inside") == 0) {
      plotoptions.inside = 1;
    }else if(strcmp(argv[i], "-nostddev") == 0
      ) {
      plotvariance = 0;
    }else if(strcmp(argv[i], "-nomod") == 0
      ) {
      plotmodindex = 0;
    }else if(strcmp(argv[i], "-f") == 0 || strcmp(argv[i], "-noflip") == 0) {
      plotoptions.doFlip = 0;
    }else if(strcmp(argv[i], "-notop") == 0) {
      plotoptions.notop = 1;
    }else if(strcmp(argv[i], "-noside") == 0) {
      plotoptions.noside = 1;
    }else if(strcmp(argv[i], "-nomain") == 0) {
      plotoptions.nomain = 1;
    }else if(strcmp(argv[i], "-overlay") == 0) {
      plotoptions.overlaypp = 1;
    }else if(strcmp(argv[i], "-noylabels") == 0) {
      plotoptions.noylabels = 1;
    }else if(strcmp(argv[i], "-noxlabels") == 0) {
      plotoptions.noxlabels = 1;
    }else if(strcmp(argv[i], "-showwedge") == 0) {
      plotoptions.showwedge = 1;
    }else if(strcmp(argv[i], "-showwedge_showmax") == 0) {
      plotoptions.showwedge_max = 1;
    }else if(strcmp(argv[i], "-normspectra") == 0) {
      plotoptions.normalise_spectra = 1;
    }else if(strcmp(argv[i], "-normside") == 0) {
      plotoptions.normaliseSide = 1;
    }else if(strcmp(argv[i], "-intnrs") == 0) {
      plotoptions.nointegrateNumbers = 0;
    }else if(strcmp(argv[i], "-intflip") == 0) {
      plotoptions.intflip = 1;
    }else if(strcmp(argv[i], "-phase") == 0) {
      plotoptions.usephase = 1;
    }else if(strcmp(argv[i], "-linestyle") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%d", &lineStyle, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-linecolor") == 0 || strcmp(argv[i], "-linecolour") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%d", &(plotoptions.lineColour), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-scalep") == 0
      ) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &profilescale, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-dlong") == 0 || strcmp(argv[i], "-dl") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &(plotoptions.dl), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-labelscale") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &(plotoptions.labelscale), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-titlescale") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &(plotoptions.titlescale), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-Imax") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &ImaxValue, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      ImaxSet = 1;
      i++;
    }else if(strcmp(argv[i], "-Imin") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &IminValue, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      IminSet = 1;
      i++;
    }else if(strcmp(argv[i], "-spmax") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &maxSubpulsePhase, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      maxSubpulsePhaseSet = 1;
      i++;
    }else if(strcmp(argv[i], "-spmin") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &minSubpulsePhase, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      minSubpulsePhaseSet = 1;
      i++;
    }else if(strcmp(argv[i], "-header") == 0) {
      i++;
    }else if(strcmp(argv[i], "-altProf") == 0) {
      altProf = i+1;
      i++;
    }else if(strcmp(argv[i], "-p3") == 0) {
      plotoptions.SelectP3Region = 1;
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &P3RegionLow, &P3RegionHigh, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-p2") == 0) {
      if(SelectP2Region == 0) {
 SelectP2Region = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &(twodfs_allinfo.f2_min), &(twodfs_allinfo.f2_max), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 twodfs2_allinfo.f2_min = twodfs_allinfo.f2_min;
 twodfs2_allinfo.f2_max = twodfs_allinfo.f2_max;
      }else {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &(twodfs2_allinfo.f2_min), &(twodfs2_allinfo.f2_max), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
      }
      i++;
    }else if(strcmp(argv[i], "-long") == 0 || strcmp(argv[i], "-l") == 0) {
      plotoptions.SelectLRegion = 1;
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &(plotoptions.LRegionLow), &(plotoptions.LRegionHigh), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-phaseslope") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f %f", &phase_slope_g, &phase_slope_o, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      do_phase_slope = 1;
      i++;
    }else if(strcmp(argv[i], "-scalel") == 0
      ) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &(plotoptions.oversaturizel), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-scale2") == 0
      ) {
      if(NrSelectedOverSaturize == 0) {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &(twodfs_allinfo.oversaturize), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 twodfs2_allinfo.oversaturize = twodfs_allinfo.oversaturize;
 NrSelectedOverSaturize = 1;
      }else {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf", &(twodfs2_allinfo.oversaturize), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
      }
      i++;
    }else if(strcmp(argv[i], "-modsigma") == 0
      ) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0,-1, "%f", &maxsigma_mod, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-stddevsigma") == 0
      ) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &maxsigma_stddev, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-modmax") == 0
      ) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &maxvalue_mod, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-stddevmax") == 0
      ) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%f", &maxvalue_stddev, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-s") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%d", &(plotoptions.ExtraVerticalMaximaSkip), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-title") == 0
      ) {
      title_index = i+1;
      i++;
    }else if(strcmp(argv[i], "-textside") == 0
      ) {
      plotoptions.textside = argv[i+1];
      i++;
    }else if(strcmp(argv[i], "-2dfsnr") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%d", &file_number, NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else if(strcmp(argv[i], "-int") == 0
      ) {
      if(twodfs_allinfo.SelectP3Integrate == 0) {
 twodfs_allinfo.SelectP3Integrate = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &(twodfs_allinfo.P3IntegrateLow), &(twodfs_allinfo.P3IntegrateHigh), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 twodfs2_allinfo.SelectP3Integrate = 1;
 twodfs2_allinfo.P3IntegrateLow = twodfs_allinfo.P3IntegrateLow;
 twodfs2_allinfo.P3IntegrateHigh = twodfs_allinfo.P3IntegrateHigh;
      }else {
 twodfs2_allinfo.SelectP3Integrate = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &(twodfs2_allinfo.P3IntegrateLow), &(twodfs2_allinfo.P3IntegrateHigh), NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
      }
      i++;
    }else if(strcmp(argv[i], "-ytop") == 0
      ) {
      plotoptions.plot_ylabeltop = 1;
    }else if(strcmp(argv[i], "-scalefig") == 0) {
      if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%lf %lf", &(plotoptions.scaleFig_x), &(plotoptions.scaleFig_y), NULL) == 0) {
 printerror(application.verbose_state.debug, "ERROR pspecFig: Cannot parse '%s' option.", argv[i]);
 return 0;
      }
      i++;
    }else {
      printerror(application.verbose_state.debug, "Unknown option: %s\nRun pspecFig without command line options for a help", argv[i]);
      terminateApplication(&application);
      return 0;
    }
  }
  copyVerboseState(application.verbose_state, &noverbose);
  noverbose.verbose = 0;
  cleanPSRData(&(twodfs_allinfo.datafile), application.verbose_state);
  cleanPSRData(&(twodfs2_allinfo.datafile), application.verbose_state);
  cleanPSRData(&lrfs, application.verbose_state);
  cleanPSRData(&AverageProfile, application.verbose_state);
  cleanPSRData(&VarianceProfile, application.verbose_state);
  cleanPSRData(&VarianceProfileErr, application.verbose_state);
  cleanPSRData(&ModProfile, application.verbose_state);
  cleanPSRData(&ModProfileErr, application.verbose_state);
  cleanPSRData(&subpulseTrackProfile, application.verbose_state);
  cleanPSRData(&subpulseTrackProfileErr, application.verbose_state);
  cleanPSRData(&subpulseAmpProfile, application.verbose_state);
  if(type_of_plots == 0) {
    if(loadLRFS(&lrfs, extprefix, longsnap, argc, argv, &plotoptions, application.verbose_state) == 0) {
      return 0;
    }
  }
  double period;
  if(loadHeaderPulseStack(&AverageProfile, lrfs, &period, type_of_plots, argc, argv, application.verbose_state) == 0) {
    return 0;
  }
  if(extprefix == 0) {
    sprintf(txt, "profile");
  }
  if(change_filename_extension(argv[argc-1], filename, txt, 1000, application.verbose_state) == 0)
    return 0;
  if(altProf > 0) {
    strncpy(filename, argv[altProf], 1000);
  }
  if(set_filename_PSRData(&AverageProfile, filename, application.verbose_state) == 0) {
    fflush(stdout);
    printerror(application.verbose_state.debug, "ERROR psecFig: Setting file name failed.");
    return 0;
  }
  if(type_of_plots == 0) {
    copy_params_PSRData(AverageProfile, &VarianceProfile, application.verbose_state);
    copy_params_PSRData(AverageProfile, &ModProfile, application.verbose_state);
    copy_params_PSRData(AverageProfile, &VarianceProfileErr, application.verbose_state);
    copy_params_PSRData(AverageProfile, &ModProfileErr, application.verbose_state);
    VarianceProfile.data = malloc(AverageProfile.NrBins*sizeof(float));
    ModProfile.data = malloc(AverageProfile.NrBins*sizeof(float));
    VarianceProfileErr.data = malloc(AverageProfile.NrBins*sizeof(float));
    ModProfileErr.data = malloc(AverageProfile.NrBins*sizeof(float));
    if(VarianceProfile.data == NULL || ModProfile.data == NULL || VarianceProfileErr.data == NULL || ModProfileErr.data == NULL) {
      printerror(application.verbose_state.debug, "Memory allocation error");
      return 0;
    }
  }else {
    copy_params_PSRData(AverageProfile, &subpulseTrackProfile, application.verbose_state);
    copy_params_PSRData(AverageProfile, &subpulseTrackProfileErr, application.verbose_state);
    copy_params_PSRData(AverageProfile, &subpulseAmpProfile, application.verbose_state);
    subpulseTrackProfile.data = malloc(AverageProfile.NrBins*sizeof(float));
    subpulseTrackProfileErr.data = malloc(AverageProfile.NrBins*sizeof(float));
    subpulseAmpProfile.data = malloc(AverageProfile.NrBins*sizeof(float));
    if(subpulseTrackProfile.data == NULL || subpulseTrackProfileErr.data == NULL || subpulseAmpProfile.data == NULL) {
      printerror(application.verbose_state.debug, "Memory allocation error");
      return 0;
    }
  }
  if(application.verbose_state.verbose)
    printf("Reading %s\n", filename);
  fout_ascii = fopen(filename, "r");
  if(fout_ascii == NULL) {
    printerror(application.verbose_state.debug, "ERROR pspecFig: Unable to open %s.", filename);
    return 0;
  }
  for(i = 0; i < AverageProfile.NrBins; i++) {
    if(type_of_plots == 0) {
      j = fscanf(fout_ascii, "%ld %f %f %f %f %f", &k, &AverageProfile.data[i], &VarianceProfile.data[i], &VarianceProfileErr.data[i], &ModProfile.data[i], &ModProfileErr.data[i]);
    }else {
      float junk;
      j = fscanf(fout_ascii, "%ld %f %f %f %f %f", &k, &AverageProfile.data[i], &junk, &junk, &junk, &junk);
    }
    if(j != 6) {
 printwarning(application.verbose_state.debug, "WARNING: It looks like profile data is rebinned? Check the units.");
 AverageProfile.fixedtsamp *= AverageProfile.NrBins/(double)(i);
 AverageProfile.tsampMode = TSAMPMODE_FIXEDTSAMP;
 AverageProfile.NrBins = i;
 printwarning(application.verbose_state.debug, "WARNING: Assuming the number of bins = %ld and the sampling time = %lf s.", AverageProfile.NrBins, AverageProfile.fixedtsamp);
 if(type_of_plots == 0) {
   copy_params_PSRData(AverageProfile, &VarianceProfile, application.verbose_state);
   copy_params_PSRData(AverageProfile, &ModProfile, application.verbose_state);
   copy_params_PSRData(AverageProfile, &VarianceProfileErr, application.verbose_state);
   copy_params_PSRData(AverageProfile, &ModProfileErr, application.verbose_state);
 }else {
   copy_params_PSRData(AverageProfile, &subpulseTrackProfile, application.verbose_state);
   copy_params_PSRData(AverageProfile, &subpulseTrackProfileErr, application.verbose_state);
   copy_params_PSRData(AverageProfile, &subpulseAmpProfile, application.verbose_state);
 }
 break;
    }
    if(k != i) {
      printerror(application.verbose_state.debug, "Unexpected bin number");
      return 0;
    }
  }
  fclose(fout_ascii);
  if(application.verbose_state.verbose)
    printf("%ld points read from %s\n", AverageProfile.NrBins, filename);
  ret = get_period(AverageProfile, 0, &period, application.verbose_state);
  if(ret == 2) {
    printerror(application.verbose_state.debug, "ERROR pspecFig (%s): Cannot obtain period", AverageProfile.filename);
    return 0;
  }
  plotoptions.fl_min = 0 + plotoptions.dl;
  plotoptions.fl_max = 360*(AverageProfile.NrBins-1)*get_tsamp(AverageProfile, 0, application.verbose_state)/period + plotoptions.dl;
  if(plotoptions.usephase) {
    plotoptions.fl_min /= 360.0;
    plotoptions.fl_max /= 360.0;
  }
  if(type_of_plots == 1) {
    if(extprefix == 0) {
      sprintf(txt, "amplitude");
    }
    if(change_filename_extension(argv[argc-1], filename, txt, 1000, application.verbose_state) == 0)
      return 0;
    if(application.verbose_state.verbose)
      printf("Reading %s\n", filename);
    fout_ascii = fopen(filename, "r");
    if(fout_ascii == NULL) {
      printerror(application.verbose_state.debug, "ERROR pspecFig: Unable to open %s.", filename);
      return 0;
    }
    for(i = 0; i < AverageProfile.NrBins; i++) {
      j = fscanf(fout_ascii, "%ld %f", &k, &subpulseAmpProfile.data[i]);
      if(j != 2) {
 printerror(application.verbose_state.debug, "Unexpected end of file (bin %d). Resolution in subpulse phase track doesn't match resolution in profile file?", i+1);
 return 0;
      }
      if(k != i) {
 printerror(application.verbose_state.debug, "Unexpected bin number");
 return 0;
      }
    }
    fclose(fout_ascii);
    if(application.verbose_state.verbose)
      printf("%ld points read from %s\n", AverageProfile.NrBins, filename);
  }
  if(type_of_plots == 1) {
    if(extprefix == 0) {
      sprintf(txt, "track");
    }
    if(change_filename_extension(argv[argc-1], filename, txt, 1000, application.verbose_state) == 0)
      return 0;
    if(application.verbose_state.verbose)
      printf("Reading %s\n", filename);
    fout_ascii = fopen(filename, "r");
    if(fout_ascii == NULL) {
      printerror(application.verbose_state.debug, "ERROR pspecFig: Unable to open %s.", filename);
      return 0;
    }
    for(i = 0; i < AverageProfile.NrBins; i++) {
      j = fscanf(fout_ascii, "%ld %f %f", &k, &subpulseTrackProfile.data[i], &subpulseTrackProfileErr.data[i]);
      if(j != 3) {
 printerror(application.verbose_state.debug, "Unexpected end of file (bin %d). Resolution in subpulse phase track doesn't match resolution in profile file?", i+1);
 return 0;
      }
      if(k != i) {
 printerror(application.verbose_state.debug, "Unexpected bin number");
 return 0;
      }
    }
    fclose(fout_ascii);
    if(application.verbose_state.verbose)
      printf("%ld points read from %s\n", AverageProfile.NrBins, filename);
  }
  twodfs_allinfo.f3_min = 0;
  twodfs_allinfo.f3_max = 0.5;
  if(type_of_plots == 0 && Load2dfs != 0) {
    if(load2dfs(&twodfs_allinfo, AverageProfile, file_number, extprefix, 0, argc, argv, plotoptions, application.verbose_state) == 0) {
      return 0;
    }
    if(LoadTwo != 0) {
      file_number++;
      if(load2dfs(&twodfs2_allinfo, AverageProfile, file_number, extprefix, 0, argc, argv, plotoptions, application.verbose_state) == 0) {
 return 0;
      }
    }
    if(plotoptions.SelectP3Region != 0) {
      twodfs_allinfo.f3_min = P3RegionLow;
      twodfs_allinfo.f3_max = P3RegionHigh;
      twodfs2_allinfo.f3_min = P3RegionLow;
      twodfs2_allinfo.f3_max = P3RegionHigh;
    }
    if(twodfs_allinfo.SelectP3Integrate == 0) {
      twodfs_allinfo.P3IntegrateLow = twodfs_allinfo.f3_min;
      twodfs_allinfo.P3IntegrateHigh = twodfs_allinfo.f3_max;
      twodfs2_allinfo.P3IntegrateLow = twodfs_allinfo.f3_min;
      twodfs2_allinfo.P3IntegrateHigh = twodfs_allinfo.f3_max;
    }
  }
  ppgopen(PlotDevice);
  ppgask(0);
  ppgslw(1);
  ppgpage();
  ppgslw(1);
  ppgscf(1);
  ppgsch(0.38*plotoptions.labelscale);
  if(plotoptions.notop == 0) {
    ppgsch(0.55*plotoptions.labelscale*plotoptions.titlescale);
    ppgscf(2);
    ppgslw(2);
    ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.15*plotoptions.scaleFig_y, 0.95);
    if(title_index > 0) {
      ppgmtxt("t",1,0.5,0.5,argv[title_index]);
    }
    ppgslw(1);
    ppgscf(1);
    ppgsch(0.38*plotoptions.labelscale);
    plotoptions.Imin = 0;
    plotoptions.Imax = 0;
    ret = get_period(AverageProfile, 0, &period, application.verbose_state);
    if(ret == 2) {
      printerror(application.verbose_state.debug, "ERROR pspecFig (%s): Cannot obtain period", AverageProfile.filename);
      return 0;
    }
    for(xi=0; xi < AverageProfile.NrBins; xi++) {
      double xpos;
      xpos = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
      if(plotoptions.usephase)
 xpos /= 360.0;
      if(xpos >= plotoptions.fl_min && xpos <= plotoptions.fl_max) {
 I = AverageProfile.data[xi];
 if(I > plotoptions.Imax)
   plotoptions.Imax = I;
 if(I < plotoptions.Imin)
   plotoptions.Imin = I;
      }
    }
    if(plotoptions.SelectLRegion != 0) {
      plotoptions.fl_min = plotoptions.LRegionLow;
      plotoptions.fl_max = plotoptions.LRegionHigh;
    }
    if(type_of_plots == 0) {
      for(xi=0; xi < ModProfile.NrBins; xi++) {
 double xpos;
 xpos = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
 if(plotoptions.usephase)
   xpos /= 360.0;
 if(xpos+plotoptions.dl >= plotoptions.fl_min && xpos+plotoptions.dl <= plotoptions.fl_max) {
   ok_flag = 1;
   if(maxsigma_mod > 0 && ModProfile.data[xi]/ModProfileErr.data[xi] < maxsigma_mod)
     ok_flag = 0;
   if(maxvalue_mod > 0 && ModProfile.data[xi] > maxvalue_mod)
     ok_flag = 0;
   if(ok_flag) {
     I = ModProfile.data[xi];
     if(I+ModProfileErr.data[xi] > plotoptions.Imax)
       plotoptions.Imax = I+ModProfileErr.data[xi];
     if(I-ModProfileErr.data[xi] < plotoptions.Imin)
       plotoptions.Imin = I-ModProfileErr.data[xi];
   }
 }
      }
      for(xi=0; xi < VarianceProfile.NrBins; xi++) {
 double xpos;
 xpos = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
 if(plotoptions.usephase)
   xpos /= 360.0;
 if(xpos + plotoptions.dl >= plotoptions.fl_min && xpos + plotoptions.dl <= plotoptions.fl_max) {
   I = VarianceProfile.data[xi];
   if(I > plotoptions.Imax)
     plotoptions.Imax = I;
   if(I < plotoptions.Imin)
     plotoptions.Imin = I;
 }
      }
    }else {
      for(xi=0; xi < subpulseAmpProfile.NrBins; xi++) {
 double xpos;
 xpos = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
 if(plotoptions.usephase)
   xpos /= 360.0;
 if(xpos + plotoptions.dl >= plotoptions.fl_min && xpos + plotoptions.dl <= plotoptions.fl_max) {
   I = subpulseAmpProfile.data[xi];
   if(I > plotoptions.Imax)
     plotoptions.Imax = I;
   if(I < plotoptions.Imin)
     plotoptions.Imin = I;
 }
      }
    }
    if(-0.05*plotoptions.Imax < plotoptions.Imin)
      plotoptions.Imin = -0.05*plotoptions.Imax;
    if(ImaxSet)
      plotoptions.Imax = ImaxValue/1.05;
    if(IminSet)
      plotoptions.Imin = IminValue/1.05;
    ppgswin(plotoptions.fl_min, plotoptions.fl_max, plotoptions.Imin, 1.05*plotoptions.Imax);
    if(LoadLRFS == 0) {
      if(plotoptions.inside) {
 if(plotoptions.noylabels) {
   ppgbox("bcnst",0.0,0,"bcst",0.0,0);
 }else {
   ppgbox("bcnst",0.0,0,"bcnst",0.0,0);
 }
      }else {
 if(plotoptions.noylabels) {
   ppgbox("bcnst",0.0,0,"bcsti",0.0,0);
 }else {
   ppgbox("bcnst",0.0,0,"bcnsti",0.0,0);
 }
      }
    }else {
      if(plotoptions.inside) {
 if(plotoptions.noylabels) {
   ppgbox("cst",0.0,0,"bcst",0.0,0);
 }else {
   ppgbox("cst",0.0,0,"bcnst",0.0,0);
 }
      }else {
 if(plotoptions.noylabels) {
   ppgbox("cst",0.0,0,"bcsti",0.0,0);
 }else {
   ppgbox("cst",0.0,0,"bcnsti",0.0,0);
 }
      }
    }
    if(plotoptions.plot_xlabel != 0 && LoadLRFS == 0) {
      ppgsch(0.3*plotoptions.labelscale);
      if(plotoptions.usephase)
 ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (phase)");
      else
 ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (deg)");
      ppgsch(0.38*plotoptions.labelscale);
    }
    ppgsls(lineStyle);
    ppgslw(2);
    ppgsci(plotoptions.lineColour);
    i = 0;
    double deg_per_sample;
    deg_per_sample = get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
    for(xi=0; xi < AverageProfile.NrBins; xi++) {
      double xpos;
      xpos = xi*deg_per_sample;
      if(plotoptions.usephase)
 xpos /= 360.0;
      if(xpos+plotoptions.dl >= plotoptions.fl_min-deg_per_sample && xpos+plotoptions.dl <= plotoptions.fl_max+deg_per_sample) {
 I = AverageProfile.data[xi]*profilescale;
 if(i == 0) {
   ppgmove(xpos+plotoptions.dl,I);
   i = 1;
 }else {
   ppgdraw(xpos+plotoptions.dl,I);
 }
      }
    }
    ppgsci(1);
    if(type_of_plots == 0) {
      if(plotvariance) {
 i = 0;
 ppgsls(1);
 ppgslw(1);
 ppgscr(20, 0.5, 0.5, 0.5);
 ppgsci(20);
 for(xi=0; xi < AverageProfile.NrBins; xi++) {
   x = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
   if(plotoptions.usephase)
     x /= 360.0;
   x += plotoptions.dl;
   if(x >= plotoptions.fl_min-deg_per_sample && x <= plotoptions.fl_max+deg_per_sample) {
     I = VarianceProfile.data[xi];
     ok_flag = 1;
     if(maxsigma_stddev > 0 && VarianceProfile.data[xi]/VarianceProfileErr.data[xi] < maxsigma_stddev)
       ok_flag = 0;
     if(maxvalue_stddev > 0 && VarianceProfile.data[xi] > maxvalue_stddev)
       ok_flag = 0;
     if(ok_flag == 0) {
       I = 0;
       i = 0;
     }else {
       if(i == 0) {
  ppgmove(x,I);
  i = 1;
       }else {
  ppgsci(20);
  ppgdraw(x,I);
       }
       ppgsci(1);
       ppgpt1(x, I, 4);
     }
   }
 }
 ppgsci(1);
      }
      if(plotmodindex) {
 i = 0;
 ppgsls(1);
 ppgslw(1);
 for(xi=0; xi < ModProfile.NrBins; xi++) {
   x = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
   if(plotoptions.usephase)
     x /= 360.0;
   x += plotoptions.dl;
   if(x >= plotoptions.fl_min-deg_per_sample && x <= plotoptions.fl_max+deg_per_sample) {
     I = ModProfile.data[xi];
     ok_flag = 1;
     if(maxsigma_mod > 0 && ModProfile.data[xi]/ModProfileErr.data[xi] < maxsigma_mod)
       ok_flag = 0;
     if(maxvalue_mod > 0 && ModProfile.data[xi] > maxvalue_mod)
       ok_flag = 0;
     if(ok_flag == 0) {
       I = 0;
       i = 0;
     }else {
       if(i == 0) {
  ppgmove(x,I);
  i = 1;
       }else {
  ppgdraw(x,I);
       }
     }
   }
 }
 ppgsls(1);
 ppgslw(1);
 for(xi=0; xi < ModProfile.NrBins; xi++) {
   x = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
   if(plotoptions.usephase)
     x /= 360.0;
   x += plotoptions.dl;
   if(x >= plotoptions.fl_min-deg_per_sample && x <= plotoptions.fl_max+deg_per_sample) {
     I = ModProfile.data[xi];
     ok_flag = 1;
     if(maxsigma_mod > 0 && ModProfile.data[xi]/ModProfileErr.data[xi] < maxsigma_mod)
       ok_flag = 0;
     if(maxvalue_mod > 0 && ModProfile.data[xi] > maxvalue_mod)
       ok_flag = 0;
     if(ok_flag != 0) {
       ppgpt1(x, I, -1);
       ppgerr1(6,x,I,ModProfileErr.data[xi], 1);
     }
   }
 }
      }
    }else {
      i = 0;
      ppgsls(4);
      ppgslw(1);
      for(xi=0; xi < AverageProfile.NrBins; xi++) {
 x = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
 if(plotoptions.usephase)
   x /= 360.0;
 x += plotoptions.dl;
 if(x >= plotoptions.fl_min-deg_per_sample && x <= plotoptions.fl_max+deg_per_sample) {
   I = subpulseAmpProfile.data[xi];
   if(i == 0) {
     ppgmove(x, I);
     i = 1;
   }else {
     ppgdraw(x, I);
   }
 }
      }
      ppgsls(1);
    }
    if(plotoptions.plot_ylabel != 0) {
      ppgsch(0.3*plotoptions.labelscale);
      if(plotoptions.plot_ylabeltop) {
 if(type_of_plots == 0) {
   ppgmtxt("l",2.8,0.5,0.5,"Intensity/Modulation index");
 }else {
   ppgmtxt("l",2.8,0.5,0.5,"Intensity");
 }
      }
      ppgsch(0.38*plotoptions.labelscale);
    }
  }
  if(type_of_plots == 0) {
    if(Load2dfs != 0 && LoadLRFS == 0) {
      char *title;
      title = NULL;
      if(title_index > 0 && plotoptions.notop && LoadLRFS == 0) {
 title = argv[title_index];
      }
      Plot2dfs(twodfs_allinfo, twodfs2_allinfo, AverageProfile, 1, 0, title, plotoptions, application.verbose_state);
      if (LoadTwo != 0) {
 Plot2dfs(twodfs_allinfo, twodfs2_allinfo, AverageProfile, 1, 1, NULL, plotoptions, application.verbose_state);
      }
    }else {
      if(LoadLRFS != 0) {
 char *title;
 title = NULL;
 if(title_index > 0 && plotoptions.notop) {
   title = argv[title_index];
 }
 PlotLRFS(lrfs, AverageProfile, twodfs_allinfo, title, plotoptions, application.verbose_state);
      }
      if(Load2dfs != 0) {
 Plot2dfs(twodfs_allinfo, twodfs2_allinfo, AverageProfile, 0, 0, NULL, plotoptions, application.verbose_state);
      }
      if(LoadTwo != 0) {
 Plot2dfs(twodfs_allinfo, twodfs2_allinfo, AverageProfile, 0, 1, NULL, plotoptions, application.verbose_state);
      }
    }
  }else {
    ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.3*plotoptions.scaleFig_y, 0.95-0.15*plotoptions.scaleFig_y);
    ppgslw(1);
    ppgscf(1);
    ppgsch(0.38*plotoptions.labelscale);
    plotoptions.Imin = 0;
    plotoptions.Imax = 0;
    for(xi=0; xi < subpulseTrackProfile.NrBins; xi++) {
      double xpos;
      xpos = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
      if(plotoptions.usephase)
 xpos /= 360.0;
      if(xpos+plotoptions.dl >= plotoptions.fl_min && xpos+plotoptions.dl <= plotoptions.fl_max) {
 I = subpulseTrackProfile.data[xi];
 if(plotoptions.doFlip)
   I *= -1;
 if(I+subpulseTrackProfileErr.data[xi] > plotoptions.Imax)
   plotoptions.Imax = I+subpulseTrackProfileErr.data[xi];
 if(I-subpulseTrackProfileErr.data[xi] < plotoptions.Imin)
   plotoptions.Imin = I-subpulseTrackProfileErr.data[xi];
      }
    }
    if(plotoptions.SelectLRegion != 0) {
      plotoptions.fl_min = plotoptions.LRegionLow;
      plotoptions.fl_max = plotoptions.LRegionHigh;
    }
    if(maxSubpulsePhaseSet)
      plotoptions.Imax = maxSubpulsePhase/1.05;
    if(minSubpulsePhaseSet)
      plotoptions.Imin = minSubpulsePhase/1.05;
    ppgswin(plotoptions.fl_min, plotoptions.fl_max, plotoptions.Imin, 1.05*plotoptions.Imax);
    if(plotoptions.inside) {
      if(plotoptions.noylabels) {
 ppgbox("bcnst",0.0,0,"bcst",0.0,0);
      }else {
 ppgbox("bcnst",0.0,0,"bcnst",0.0,0);
      }
    }else {
      if(plotoptions.noylabels) {
 ppgbox("bcnst",0.0,0,"bcsti",0.0,0);
      }else {
 ppgbox("bcnst",0.0,0,"bcnsti",0.0,0);
      }
    }
    if(plotoptions.plot_xlabel != 0 && LoadLRFS == 0) {
      ppgsch(0.3*plotoptions.labelscale);
      if(plotoptions.usephase)
 ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (phase)");
      else
 ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (deg)");
      ppgsch(0.38*plotoptions.labelscale);
    }
    if(do_phase_slope) {
      ppgsls(4);
      if(plotoptions.usephase == 0) {
 x = plotoptions.dl;
 I = x*phase_slope_g + phase_slope_o;
 I = derotate_deg(I);
 ppgmove(x, I);
 ppgdraw(x+360, I+phase_slope_g*360);
 ppgmove(x, I-360);
 ppgdraw(x+360, I-360+phase_slope_g*360);
 ppgmove(x, I-2*360);
 ppgdraw(x+360, I-2*360+phase_slope_g*360);
 ppgmove(x, I+360);
 ppgdraw(x+360, I+360+phase_slope_g*360);
 ppgmove(x, I+2*360);
 ppgdraw(x+360, I+2*360+phase_slope_g*360);
      }else {
 x = plotoptions.dl;
 I = x*360.0*phase_slope_g + phase_slope_o;
 I = derotate_deg(I);
 ppgmove(x, I);
 ppgdraw(x+1.0, I+phase_slope_g*360);
 ppgmove(x, I-360);
 ppgdraw(x+1.0, I-360+phase_slope_g*360);
 ppgmove(x, I-2*360);
 ppgdraw(x+1.0, I-2*360+phase_slope_g*360);
 ppgmove(x, I+360);
 ppgdraw(x+1.0, I+360+phase_slope_g*360);
 ppgmove(x, I+2*360);
 ppgdraw(x+1.0, I+2*360+phase_slope_g*360);
      }
      ppgsls(1);
    }
    ppgscr(20, 0.5, 0.5, 0.5);
    ppgsci(20);
    ppgsls(1);
    ppgslw(1);
    for(xi=0; xi < subpulseTrackProfile.NrBins; xi++) {
      x = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
      if(plotoptions.usephase)
 x /= 360.0;
      x += plotoptions.dl;
      if(x >= plotoptions.fl_min && x <= plotoptions.fl_max) {
 I = subpulseTrackProfile.data[xi];
 if(plotoptions.doFlip)
   I *= -1;
 ppgerr1(6,x,I,subpulseTrackProfileErr.data[xi], 1);
 ppgerr1(6,x,I+360,subpulseTrackProfileErr.data[xi], 1);
 ppgerr1(6,x,I-360,subpulseTrackProfileErr.data[xi], 1);
      }
    }
    ppgsci(1);
    ppgslw(3);
    for(xi=0; xi < subpulseTrackProfile.NrBins; xi++) {
      x = xi*get_tsamp(AverageProfile, 0, application.verbose_state)*360.0/period;
      if(plotoptions.usephase)
 x /= 360.0;
      x += plotoptions.dl;
      if(x >= plotoptions.fl_min && x <= plotoptions.fl_max) {
 I = subpulseTrackProfile.data[xi];
 if(plotoptions.doFlip)
   I *= -1;
 ppgpt1(x, I, -1);
 ppgpt1(x, I+360, -1);
 ppgpt1(x, I-360, -1);
      }
    }
    ppgslw(1);
    if(plotoptions.plot_xlabel != 0) {
      ppgsch(0.3*plotoptions.labelscale);
      if(plotoptions.usephase)
 ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (phase)");
      else
 ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (deg)");
      ppgsch(0.38*plotoptions.labelscale);
    }
    if(plotoptions.plot_ylabel != 0) {
      ppgsch(0.3*plotoptions.labelscale);
      if(plotoptions.usephase)
 ppgmtxt("l",2.8,0.5,0.5,"Subpulse phase (phase)");
      else
 ppgmtxt("l",2.8,0.5,0.5,"Subpulse phase (deg)");
      ppgsch(0.38*plotoptions.labelscale);
    }
  }
  if(preprocess_checknan(AverageProfile, 1, noverbose)) {
    printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
  }
  if(preprocess_checkinf(AverageProfile, 1, noverbose)) {
    printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
  }
  if(type_of_plots == 0) {
    if(preprocess_checknan(VarianceProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(VarianceProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checknan(ModProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(ModProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checknan(VarianceProfileErr, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(VarianceProfileErr, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checknan(ModProfileErr, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(ModProfileErr, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
  }else {
    if(preprocess_checknan(subpulseTrackProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(subpulseTrackProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checknan(subpulseTrackProfileErr, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(subpulseTrackProfileErr, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checknan(subpulseAmpProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have NAN's in them. Artifacts can be expected in the plot.");
    }
    if(preprocess_checkinf(subpulseAmpProfile, 1, noverbose)) {
      printwarning(application.verbose_state.debug, "WARNING: The profile data appears to have INF's in them. Artifacts can be expected in the plot.");
    }
  }
  ppgend();
  closePSRData(&AverageProfile, 0, application.verbose_state);
  closePSRData(&(twodfs_allinfo.datafile), 0, application.verbose_state);
    closePSRData(&(twodfs2_allinfo.datafile), 0, application.verbose_state);
  closePSRData(&lrfs, 0, application.verbose_state);
  closePSRData(&VarianceProfile, 0, application.verbose_state);
  closePSRData(&ModProfile, 0, application.verbose_state);
  closePSRData(&VarianceProfileErr, 0, application.verbose_state);
  closePSRData(&ModProfileErr, 0, application.verbose_state);
  closePSRData(&subpulseTrackProfile, 0, application.verbose_state);
  closePSRData(&subpulseTrackProfileErr, 0, application.verbose_state);
  closePSRData(&subpulseAmpProfile, 0, application.verbose_state);
  terminateApplication(&application);
  return 0;
}
void GetExtremesSubsetVertical(twodfs_def twodfs_allinfo, float *Imin, float *Imax)
{
  float I, x, y;
  int xi, yi;
  *Imin = 0;
  *Imax = 0;
  for(xi = 0; xi < twodfs_allinfo.datafile.NrBins; xi++) {
    I = 0;
    for(yi = 0; yi < twodfs_allinfo.datafile.NrSubints; yi++) {
      pgplotMapCoordinateInverse(&x, &y, xi, yi);
      if(y >= twodfs_allinfo.P3IntegrateLow && y <= twodfs_allinfo.P3IntegrateHigh) {
 if(x >= twodfs_allinfo.f2_min && x <= twodfs_allinfo.f2_max) {
   I += twodfs_allinfo.datafile.data[yi*twodfs_allinfo.datafile.NrBins+xi];
 }
      }
    }
    if(I < *Imin)
      *Imin = I;
    if(I > *Imax)
      *Imax = I;
  }
  *Imin *= 2.0;
  *Imax *= 2.0;
}
void GetExtremesSubsetHorizontal(datafile_definition spectrum, int istwodfs, double f2_min, double f2_max, double f3_min, double f3_max, float *Imin, float *Imax, plotoptions_def plotoptions)
{
  float I, x, y;
  int xi, yi, ok;
  *Imin = 0;
  *Imax = 0;
  for(yi = plotoptions.ExtraVerticalMaximaSkip; yi < spectrum.NrSubints; yi++) {
    I = 0;
    for(xi = 0; xi < spectrum.NrBins; xi++) {
      pgplotMapCoordinateInverse(&x, &y, xi, yi);
      ok = 0;
      if(istwodfs == 0) {
 if(x >= plotoptions.fl_min && x <= plotoptions.fl_max && y >= f3_min && y <= f3_max) {
   ok = 1;
 }
      }else {
 if(x >= f2_min && x <= f2_max && y >= f3_min && y <= f3_max) {
   ok = 1;
 }
      }
      if(ok) {
 I += spectrum.data[yi*spectrum.NrBins+xi];
      }
    }
    if(I < *Imin)
      *Imin = I;
    if(I > *Imax)
      *Imax = I;
  }
  *Imin *= 2.0;
  *Imax *= 2.0;
}
void IntegrateSubsetHorizontal(datafile_definition lrfs, twodfs_def twodfs_allinfo, twodfs_def twodfs2_allinfo, int twodfsonly, int Number, int normalise, plotoptions_def plotoptions)
{
  float I, Imin, Imax, x, y;
  int xi, yi;
  float offset;
  offset = 0+0.03*(plotoptions.labelscale-1.0)*((float)Number+1.0);
  if(Number == -1) {
    ppgsvp(0.2-0.04*plotoptions.scaleFig_x, 0.2, 0.95-0.3*plotoptions.scaleFig_y, 0.95-0.15*plotoptions.scaleFig_y);
  }else {
    if(twodfsonly == 0) {
      ppgsvp(0.2-0.04*plotoptions.scaleFig_x, 0.2, 0.95-0.48*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y, 0.95-0.33*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y);
    }else {
      ppgsvp(0.2-0.04*plotoptions.scaleFig_x, 0.2, 0.95-0.33*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y, 0.95-0.18*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y);
    }
  }
  if(Number == -1) {
    GetExtremesSubsetHorizontal(lrfs, 0, 0.0, 0.0, twodfs_allinfo.f3_min, twodfs_allinfo.f3_max, &Imin, &Imax, plotoptions);
  }else if(Number == 0) {
    GetExtremesSubsetHorizontal(twodfs_allinfo.datafile, 1, twodfs_allinfo.f2_min, twodfs_allinfo.f2_max, twodfs_allinfo.f3_min, twodfs_allinfo.f3_max, &Imin, &Imax, plotoptions);
  }else {
    GetExtremesSubsetHorizontal(twodfs2_allinfo.datafile, 1, twodfs2_allinfo.f2_min, twodfs2_allinfo.f2_max, twodfs2_allinfo.f3_min, twodfs2_allinfo.f3_max, &Imin, &Imax, plotoptions);
  }
  double scale;
  scale = fabs(Imax);
  if(fabs(Imin) > scale)
    scale = fabs(Imin);
  if(normalise == 0)
    scale = 1.0;
  ppgswin(Imin/scale,1.05*Imax/scale, twodfs_allinfo.f3_min,twodfs_allinfo.f3_max);
  ppgbox("bcv",0.0,0,"bnsti",0.0,0);
  ppgbox("",0.0,0,"c",0.0,0);
  if(plotoptions.textside != NULL) {
    if(twodfsonly == 0) {
      if(Number == -1) {
 ppgsch(0.55*plotoptions.labelscale*plotoptions.titlescale);
 ppgscf(2);
 ppgslw(2);
 if (plotoptions.notop == 1) {
   ppgtext(Imin, (twodfs_allinfo.f3_max-twodfs_allinfo.f3_min)*1.15, plotoptions.textside);
 } else {
   ppgtext(Imin, (twodfs_allinfo.f3_max-twodfs_allinfo.f3_min)*2.095, plotoptions.textside);
 }
 ppgslw(1);
 ppgscf(1);
 ppgsch(0.38*plotoptions.labelscale);
      }
    }else {
      if(strlen(plotoptions.textside) != 0 && plotoptions.notop == 0) {
 ppgsch(0.55*plotoptions.labelscale*plotoptions.titlescale);
 ppgscf(2);
 ppgslw(2);
 ppgtext(Imin, (twodfs_allinfo.f3_max-twodfs_allinfo.f3_min)+0.65, plotoptions.textside);
 ppgslw(1);
 ppgscf(1);
 ppgsch(0.38*plotoptions.labelscale);
      }else if (strlen(plotoptions.textside) != 0 && plotoptions.notop == 1) {
 ppgsch(0.55*plotoptions.labelscale*plotoptions.titlescale);
 ppgscf(2);
 ppgslw(2);
 ppgtext(Imin, (twodfs_allinfo.f3_max-twodfs_allinfo.f3_min)+0.1, plotoptions.textside);
 ppgslw(1);
 ppgscf(1);
 ppgsch(0.38*plotoptions.labelscale);
      }
    }
  }
  y = floor(log10(Imax*0.33/scale));
  x = floor(Imax*0.33/(pow(10,y)*scale));
  x = x*pow(10,y);
  char labelnumbers[3];
  if(plotoptions.nointegrateNumbers == 0) {
    strcpy(labelnumbers, "n");
  }else {
    strcpy(labelnumbers, "");
  }
  ppgsch(0.38*plotoptions.labelscale*0.66);
  ppgaxis(labelnumbers,0,twodfs_allinfo.f3_max,Imax*0.33/scale,twodfs_allinfo.f3_max,0,Imax*0.33/scale,x,1,0.3,0,0,-0.5,90);
  ppgsch(0.38*plotoptions.labelscale);
  datafile_definition *spectrum;
  double xmin, xmax, ymin, ymax;
  if(Number == -1) {
    spectrum = &lrfs;
    xmin = plotoptions.fl_min;
    xmax = plotoptions.fl_max;
    ymin = twodfs_allinfo.f3_min;
    ymax = twodfs_allinfo.f3_max;
  }else if(Number == 0) {
    spectrum = &(twodfs_allinfo.datafile);
    xmin = twodfs_allinfo.f2_min;
    xmax = twodfs_allinfo.f2_max;
    ymin = twodfs_allinfo.f3_min;
    ymax = twodfs_allinfo.f3_max;
  }else if(Number == 1) {
    spectrum = &(twodfs2_allinfo.datafile);
    xmin = twodfs2_allinfo.f2_min;
    xmax = twodfs2_allinfo.f2_max;
    ymin = twodfs2_allinfo.f3_min;
    ymax = twodfs2_allinfo.f3_max;
  }
  for(yi = 0; yi < spectrum->NrSubints; yi++) {
    I = 0;
    for(xi = 0; xi < spectrum->NrBins; xi++) {
      pgplotMapCoordinateInverse(&x, &y, xi, yi);
      if(x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
 I += spectrum->data[yi*spectrum->NrBins+xi];
      }
    }
    if(yi == 0) {
      ppgmove(2.0*I/scale, y);
    }else {
      ppgdraw(2.0*I/scale, y);
    }
  }
}
void PlotLRFS(datafile_definition lrfs, datafile_definition AverageProfile, twodfs_def twodfs_allinfo, char *title, plotoptions_def plotoptions, verbose_definition verbose)
{
  float I;
  int xi, i;
  ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.3*plotoptions.scaleFig_y, 0.95-0.15*plotoptions.scaleFig_y);
  if(title != NULL) {
    ppgsch(0.55*plotoptions.labelscale*plotoptions.titlescale);
    ppgscf(2);
    ppgslw(2);
    ppgmtxt("t",1,0.5,0.5,title);
    ppgslw(1);
    ppgscf(1);
    ppgsch(0.38*plotoptions.labelscale);
  }
  if(plotoptions.SelectLRegion != 0) {
    plotoptions.fl_min = plotoptions.LRegionLow;
    plotoptions.fl_max = plotoptions.LRegionHigh;
  }
  ppgswin(plotoptions.fl_min,plotoptions.fl_max,twodfs_allinfo.f3_min,twodfs_allinfo.f3_max);
  if(plotoptions.plot_xlabel != 0) {
    ppgsch(0.3*plotoptions.labelscale);
    if(plotoptions.usephase)
      ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (phase)");
    else
      ppgmtxt("b",3.0,0.5,0.5,"Pulse longitude (deg)");
    ppgsch(0.38*plotoptions.labelscale);
  }
  if(plotoptions.plot_ylabel != 0) {
    ppgsch(0.3*plotoptions.labelscale);
    if(plotoptions.noside) {
      if(plotoptions.plot_ylabel == 1)
 ppgmtxt("l",2.8,0.5,0.5,"Fluctuation frequency (cpp)");
      else if(plotoptions.plot_ylabel == 2)
 ppgmtxt("l",2.8,0.5,0.5,"Fluctuation frequency (P\\d0\\u/P\\d3\\u)");
      else if(plotoptions.plot_ylabel == 3)
 ppgmtxt("l",2.8,0.5,0.5,"Fluctuation frequency (P/P\\d3\\u)");
    }else {
      float offset;
      offset = 10 - (plotoptions.labelscale-1.0)*(7-2.8);
      offset += 7*(plotoptions.scaleFig_x-1.0);
      if(plotoptions.plot_ylabel == 1)
 ppgmtxt("l",offset,0.5,0.5,"Fluctuation frequency (cpp)");
      else if(plotoptions.plot_ylabel == 2)
 ppgmtxt("l",offset,0.5,0.5,"Fluctuation frequency (P\\d0\\u/P\\d3\\u)");
      else if(plotoptions.plot_ylabel == 3)
 ppgmtxt("l",offset,0.5,0.5,"Fluctuation frequency (P/P\\d3\\u)");
    }
    ppgsch(0.38*plotoptions.labelscale);
  }
  pgplot_options_definition pgplot_options;
  pgplot_clear_options(&pgplot_options);
  pgplot_options.box.box_labelsize = 0.3*plotoptions.labelscale;
  pgplot_options.viewport.dontopen = 1;
  pgplot_options.viewport.dontclose = 1;
  pgplot_options.viewport.noclear = 1;
  double period;
  int ret;
  ret = get_period(lrfs, 0, &period, verbose);
  if(ret == 2) {
    printerror(verbose.debug, "ERROR pspecFig (%s): Cannot obtain period", lrfs.filename);
    exit(0);
  }
  double xright;
  xright = 360*(lrfs.NrBins-1)*get_tsamp(lrfs, 0, verbose)/period;
  if(plotoptions.usephase)
    xright /= 360.0;
  pgplotMap(&pgplot_options, lrfs.data, lrfs.NrBins, lrfs.NrSubints, 0+plotoptions.dl, xright+plotoptions.dl, plotoptions.fl_min, plotoptions.fl_max, 0, 0.5, twodfs_allinfo.f3_min, twodfs_allinfo.f3_max, PPGPLOT_GRAYSCALE, 0, 0, 0, NULL, 0, 0, plotoptions.oversaturizel, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, plotoptions.showwedge, plotoptions.showwedge_max, 0, 0, verbose);
  ret = get_period(AverageProfile, 0, &period, verbose);
  if(ret == 2) {
    printerror(verbose.debug, "ERROR pspecFig (%s): Cannot obtain period", AverageProfile.filename);
    exit(0);
  }
  if(plotoptions.overlaypp) {
    ppgsci(plotoptions.lineColour);
    ppgsls(1);
    ppgslw(1);
    i = 0;
    for(xi=0; xi < AverageProfile.NrBins; xi++) {
      double xpos;
      xpos = xi*get_tsamp(AverageProfile, 0, verbose)*360.0/period;
      if(plotoptions.usephase)
 xpos /= 360.0;
      if(xpos+plotoptions.dl >= plotoptions.fl_min && xpos+plotoptions.dl <= plotoptions.fl_max) {
 I = AverageProfile.data[xi]*0.5;
 if(i == 0) {
   ppgmove(xpos+plotoptions.dl,I);
   i = 1;
 }else {
   ppgdraw(xpos+plotoptions.dl,I);
 }
      }
    }
    ppgsci(1);
  }
  ppgswin(plotoptions.fl_min,plotoptions.fl_max,twodfs_allinfo.f3_min,twodfs_allinfo.f3_max);
  ppgsch(0.38*plotoptions.labelscale);
  if(plotoptions.noside) {
    if(plotoptions.inside) {
      if(plotoptions.noylabels)
 ppgbox("bcnst",0.0,0,"bcst",0.0,0);
      else
 ppgbox("bcnst",0.0,0,"bcnst",0.0,0);
    }else {
      if(plotoptions.noylabels)
 ppgbox("bcnsti",0.0,0,"bcsti",0.0,0);
      else
 ppgbox("bcnsti",0.0,0,"bcnsti",0.0,0);
    }
  }else {
    if(plotoptions.inside) {
      ppgbox("bcnst",0.0,0,"cst",0.0,0);
    }else {
      ppgbox("bcnsti",0.0,0,"csti",0.0,0);
    }
  }
  if(plotoptions.noside == 0) {
    IntegrateSubsetHorizontal(lrfs, twodfs_allinfo, twodfs_allinfo, 0, -1, plotoptions.normaliseSide, plotoptions);
  }
  ppgsch(0.38*plotoptions.labelscale);
}
void IntegrateSubsetVertical(twodfs_def twodfs_allinfo, int twodfsonly, int Number, int normalise, int nointegrateNumbers, plotoptions_def plotoptions)
{
  float I, Imin, Imax, x, y, offset;
  int xi, yi;
  offset = 0+0.03*(plotoptions.labelscale-1.0)*((float)Number+1.0);
  if(twodfsonly == 0) {
    ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.52*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y, 0.95-0.48*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y);
  }else {
    ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.37*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y, 0.95-0.33*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y);
  }
  GetExtremesSubsetVertical(twodfs_allinfo, &Imin, &Imax);
  double scale;
  scale = fabs(Imax);
  if(fabs(Imin) > scale)
    scale = fabs(Imin);
  if(normalise == 0)
    scale = 1.0;
  ppgswin(twodfs_allinfo.f2_min,twodfs_allinfo.f2_max,Imin/scale,1.05*Imax/scale);
  if(plotoptions.noxlabels == 0)
    ppgbox("bnst",0.0,0,"bcvi",0.0,0);
  else
    ppgbox("bst",0.0,0,"bcvi",0.0,0);
  ppgbox("c",0.0,0,"",0.0,0);
  y = floor(log10(Imax*0.7/scale));
  x = floor(Imax*0.7/(pow(10,y)*scale));
  x = x*pow(10,y);
  char labelnumbers[3];
  if(plotoptions.nointegrateNumbers == 0) {
    strcpy(labelnumbers, "n");
  }else {
    strcpy(labelnumbers, "");
  }
  ppgsch(0.38*plotoptions.labelscale*0.66);
  ppgaxis(labelnumbers,twodfs_allinfo.f2_min, 0,twodfs_allinfo.f2_min, Imax*.7/scale,0,Imax*.7/scale,x,1,0.3,0,0,-0.8,90);
  ppgsch(0.38*plotoptions.labelscale);
  int direction;
  for(direction = 1-plotoptions.intflip; direction < 2; direction++) {
    if(direction == 0) {
      ppgsci(8);
      ppgsls(2);
    }else {
      ppgsci(1);
      ppgsls(1);
    }
    for(xi = 0; xi < twodfs_allinfo.datafile.NrBins; xi++) {
      I = 0;
      for(yi = 0; yi < twodfs_allinfo.datafile.NrSubints; yi++) {
 pgplotMapCoordinateInverse(&x, &y, xi, yi);
 if(y >= twodfs_allinfo.P3IntegrateLow && y <= twodfs_allinfo.P3IntegrateHigh)
   I += twodfs_allinfo.datafile.data[yi*twodfs_allinfo.datafile.NrBins+xi];
      }
      if(direction == 0)
 x *= -1.0;
      if(xi == 0) {
 ppgmove(x,2.0*I/scale);
      }else {
 ppgdraw(x,2.0*I/scale);
      }
    }
  }
}
void Plot2dfs(twodfs_def twodfs_allinfo, twodfs_def twodfs2_allinfo, datafile_definition AverageProfile, int twodfsonly, int Number, char *title, plotoptions_def plotoptions, verbose_definition verbose)
{
  float offset;
  offset = 0+0.03*(plotoptions.labelscale-1.0)*((float)Number+1.0);
  if(twodfsonly == 0) {
    ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.48*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y, 0.95-0.33*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y);
  }else {
    ppgsvp(0.2, 0.2+0.11*plotoptions.scaleFig_x, 0.95-0.33*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y, 0.95-0.18*plotoptions.scaleFig_y-0.22*Number*plotoptions.scaleFig_y-offset*plotoptions.scaleFig_y);
  }
  if(twodfsonly && title != NULL) {
    ppgsch(0.55*plotoptions.labelscale*plotoptions.titlescale);
    ppgscf(2);
    ppgslw(2);
    ppgmtxt("t",1,0.5,0.5,title);
    ppgslw(1);
    ppgscf(1);
    ppgsch(0.38*plotoptions.labelscale);
  }
  ppgswin(twodfs_allinfo.f2_min,twodfs_allinfo.f2_max,twodfs_allinfo.f3_min,twodfs_allinfo.f3_max);
  if(plotoptions.plot_ylabel != 0) {
    ppgsch(0.3*plotoptions.labelscale);
    if(plotoptions.noside) {
      if(plotoptions.plot_ylabel == 1)
 ppgmtxt("l",2.8,0.5,0.5,"Fluctuation frequency (cpp)");
      else if(plotoptions.plot_ylabel == 2)
 ppgmtxt("l",2.8,0.5,0.5,"Fluctuation frequency (P\\d0\\u/P\\d3\\u)");
      else if(plotoptions.plot_ylabel == 3)
 ppgmtxt("l",2.8,0.5,0.5,"Fluctuation frequency (P/P\\d3\\u)");
    }else {
      float offset;
      offset = 10 - (plotoptions.labelscale-1.0)*(7-2.8);
      offset += 7*(plotoptions.scaleFig_x-1.0);
      if(plotoptions.plot_ylabel == 1)
 ppgmtxt("l",offset,0.5,0.5,"Fluctuation frequency (cpp)");
      else if(plotoptions.plot_ylabel == 2)
 ppgmtxt("l",offset,0.5,0.5,"Fluctuation frequency (P\\d0\\u/P\\d3\\u)");
      else if(plotoptions.plot_ylabel == 3)
 ppgmtxt("l",offset,0.5,0.5,"Fluctuation frequency (P/P\\d3\\u)");
    }
    ppgsch(0.38*plotoptions.labelscale);
  }
  pgplot_options_definition pgplot_options;
  pgplot_clear_options(&pgplot_options);
  pgplot_options.box.box_labelsize = 0.3*plotoptions.labelscale;
  pgplot_options.viewport.noclear = 1;
  pgplot_options.viewport.dontopen = 1;
  pgplot_options.viewport.dontclose = 1;
  if(!plotoptions.nomain) {
    datafile_definition *curtwodfs;
    float cur_f2_min, cur_f2_max;
    float cur_f3_min, cur_f3_max;
    float cur_oversaturize;
    if(Number == 0) {
      curtwodfs = &(twodfs_allinfo.datafile);
      cur_f2_min = twodfs_allinfo.f2_min;
      cur_f2_max = twodfs_allinfo.f2_max;
      cur_f3_min = twodfs_allinfo.f3_min;
      cur_f3_max = twodfs_allinfo.f3_max;
      cur_oversaturize = twodfs_allinfo.oversaturize;
    }else {
      curtwodfs = &(twodfs2_allinfo.datafile);
      cur_f2_min = twodfs2_allinfo.f2_min;
      cur_f2_max = twodfs2_allinfo.f2_max;
      cur_f3_min = twodfs2_allinfo.f3_min;
      cur_f3_max = twodfs2_allinfo.f3_max;
      cur_oversaturize = twodfs2_allinfo.oversaturize;
    }
    float cpp_min;
    float cpp_max;
    if(curtwodfs->xrangeset) {
      cpp_min = curtwodfs->xrange[0];
      cpp_max = curtwodfs->xrange[1];
    }else {
      cpp_min = -AverageProfile.NrBins/2.0-0.5*AverageProfile.NrBins/(float)curtwodfs->NrBins;
      cpp_max = +AverageProfile.NrBins/2.0-0.5*AverageProfile.NrBins/(float)curtwodfs->NrBins;
    }
    if(plotoptions.doFlip) {
      cpp_min *= -1.0;
      cpp_max *= -1.0;
    }
    pgplotMap(&pgplot_options, curtwodfs->data, curtwodfs->NrBins, curtwodfs->NrSubints, cpp_min, cpp_max, cur_f2_min, cur_f2_max, 0, 0.5, cur_f3_min, cur_f3_max, PPGPLOT_GRAYSCALE, 0, 0, 0, NULL, 0, 0, cur_oversaturize, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, plotoptions.showwedge, plotoptions.showwedge_max, 0, 0, verbose);
    ppgsch(0.38*plotoptions.labelscale);
    if(plotoptions.noside) {
      if(plotoptions.inside) {
 if(plotoptions.noylabels)
   ppgbox("cst",0.0,0,"bcst",0.0,0);
 else
   ppgbox("cst",0.0,0,"bcnst",0.0,0);
      }else {
 if(plotoptions.noylabels)
   ppgbox("csti",0.0,0,"bcsti",0.0,0);
 else
   ppgbox("csti",0.0,0,"bcnsti",0.0,0);
      }
    }else {
      if(plotoptions.inside) {
 ppgbox("cst",0.0,0,"cst",0.0,0);
      }else {
 ppgbox("csti",0.0,0,"csti",0.0,0);
      }
    }
    ppgbox("c",0.0,0,"c",0.0,0);
    if(twodfs_allinfo.SelectP3Integrate != 0) {
      ppgsls(2);
      ppgslw(1);
      if(Number == 0) {
 ppgmove(twodfs_allinfo.f2_min, twodfs_allinfo.P3IntegrateLow);
 ppgdraw(twodfs_allinfo.f2_max, twodfs_allinfo.P3IntegrateLow);
      }else {
 ppgmove(twodfs_allinfo.f2_min, twodfs2_allinfo.P3IntegrateLow);
 ppgdraw(twodfs_allinfo.f2_max, twodfs2_allinfo.P3IntegrateLow);
      }
      ppgsls(2);
      if(Number == 0) {
 ppgmove(twodfs_allinfo.f2_min, twodfs_allinfo.P3IntegrateHigh);
 ppgdraw(twodfs_allinfo.f2_max, twodfs_allinfo.P3IntegrateHigh);
      }else {
 ppgmove(twodfs_allinfo.f2_min, twodfs2_allinfo.P3IntegrateHigh);
 ppgdraw(twodfs_allinfo.f2_max, twodfs2_allinfo.P3IntegrateHigh);
      }
      ppgsls(1);
      ppgslw(1);
    }
  }
  if(Number == 0) {
    IntegrateSubsetVertical(twodfs_allinfo, twodfsonly, Number, plotoptions.normaliseSide, plotoptions.nointegrateNumbers, plotoptions);
  }else {
    IntegrateSubsetVertical(twodfs2_allinfo, twodfsonly, Number, plotoptions.normaliseSide, plotoptions.nointegrateNumbers, plotoptions);
  }
  if(plotoptions.plot_xlabel == 1) {
    ppgsch(0.3*plotoptions.labelscale);
    ppgmtxt("b",3.0,0.5,0.5,"Fluctuation frequency (cpp)");
    ppgsch(0.38*plotoptions.labelscale);
  }else if(plotoptions.plot_xlabel == 2) {
    ppgsch(0.3*plotoptions.labelscale);
    ppgmtxt("b",3.0,0.5,0.5,"Fluctuation frequency (P\\d0\\u/P\\d2\\u)");
    ppgsch(0.38*plotoptions.labelscale);
  }else if(plotoptions.plot_xlabel == 3) {
    ppgsch(0.3*plotoptions.labelscale);
    ppgmtxt("b",3.0,0.5,0.5,"Fluctuation frequency (P/P\\d2\\u)");
    ppgsch(0.38*plotoptions.labelscale);
  }
  if(!(plotoptions.noside)) {
    IntegrateSubsetHorizontal(twodfs_allinfo.datafile, twodfs_allinfo, twodfs2_allinfo, twodfsonly, Number, plotoptions.normaliseSide, plotoptions);
  }
}
int loadLRFS(datafile_definition *lrfs, int extprefix, int longsnap, int argc, char **argv, plotoptions_def *plotoptions, verbose_definition verbose)
{
  char filename[MaxFilenameLength], txt[MaxFilenameLength];
  datafile_definition clone;
  if(extprefix == 0) {
    sprintf(txt, "lrfs");
  }
  if(change_filename_extension(argv[argc-1], filename, txt, MaxFilenameLength, verbose) == 0) {
    return 0;
  }
  if(verbose.verbose) {
    printf("Reading %s\n", filename);
  }
  closePSRData(lrfs, 0, verbose);
  if(!openPSRData(lrfs, filename, 0, 0, 1, 0, verbose)) {
    return 0;
  }
  if(PSRDataHeader_parse_commandline(lrfs, argc, argv, verbose) == 0) {
    return 0;
  }
  double period;
  int ret;
  ret = get_period(*lrfs, 0, &period, verbose);
  if(ret == 2) {
    printerror(verbose.debug, "ERROR pspecFig (%s): Cannot obtain period", lrfs->filename);
    return 0;
  }
  if(period < 0.001) {
    printerror(verbose.debug, "ERROR pspecFig (%s): The period does not appear to be set in the header. Consider using the -header option.", lrfs->filename);
    return 0;
  }
  if(get_tsamp(*lrfs, 0, verbose) < 0.0000001) {
    printerror(verbose.debug, "ERROR pspecFig (%s): The sampling time does not appear to be set in the header. Consider using the -header option.", lrfs->filename);
    return 0;
  }
  if(verbose.verbose) {
    printf("%ldx%ld points read from lrfs\n", lrfs->NrBins, lrfs->NrSubints);
  }
  if(lrfs->NrPols > 1) {
    if(preprocess_polselect(*lrfs, &clone, 0, verbose) == 0) {
      printerror(verbose.debug, "ERROR pspecFig (%s): Error selecting polarization channel 0.", lrfs->filename);
      return 0;
    }
    swap_orig_clone(lrfs, &clone, verbose);
  }
  long i;
  double max;
  if(plotoptions->normalise_spectra) {
    for(i = 0; i < lrfs->NrSubints * lrfs->NrBins; i++) {
      if(i == 0 || fabs(lrfs->data[i]) > max) {
 max = fabs(lrfs->data[i]);
      }
    }
    for(i = 0; i < lrfs->NrSubints * lrfs->NrBins; i++) {
      lrfs->data[i] /= 0.999*max;
    }
  }
  return 1;
}
int loadHeaderPulseStack(datafile_definition *AverageProfile, datafile_definition lrfs, double *period, int type_of_plots, int argc, char **argv, verbose_definition verbose)
{
  if(verbose.verbose)
    printf("Reading %s\n", argv[argc-1]);
  closePSRData(AverageProfile, 0, verbose);
  if(!openPSRData(AverageProfile, argv[argc-1], 0, 0, 0, 0, verbose))
    return 0;
  if(!readHeaderPSRData(AverageProfile, 0, 0, verbose))
    return 0;
  if(PSRDataHeader_parse_commandline(AverageProfile, argc, argv, verbose) == 0)
    return 0;
  if(type_of_plots == 0) {
    if(AverageProfile->NrBins != lrfs.NrBins) {
      printwarning(verbose.debug, "WARNING: Nr of bins in pulse stack and LRFS do not match. It looks like data is rebinned? Check the units.");
    }
  }
  int ret_prd;
  ret_prd = get_period(*AverageProfile, 0, period, verbose);
  if(ret_prd == 2) {
    printerror(verbose.debug, "ERROR pspecFig (%s): Cannot obtain period", AverageProfile->filename);
    return 0;
  }
  if(*period < 0.001) {
    printerror(verbose.debug, "ERROR pspecFig (%s): The period does not appear to be set in the header. Consider using the -header option.", AverageProfile->filename);
    return 0;
  }
  if(get_tsamp(*AverageProfile, 0, verbose) < 0.0000001) {
    printerror(verbose.debug, "ERROR pspecFig (%s): The sampling time does not appear to be set in the header. Consider using the -header option.", AverageProfile->filename);
    return 0;
  }
  closePSRData(AverageProfile, 1, verbose);
  AverageProfile->format = MEMORY_format;
  AverageProfile->NrSubints = 1;
  AverageProfile->NrFreqChan = 1;
  AverageProfile->NrPols = 1;
  AverageProfile->data = malloc(AverageProfile->NrBins*sizeof(float));
  if(AverageProfile->data == NULL) {
    printerror(verbose.debug, "ERROR pspecFig (%s): Memory allocation error", AverageProfile->filename);
    return 0;
  }
  return 1;
}
int load2dfs(twodfs_def *twodfs_allinfo, datafile_definition AverageProfile, int file_number, int extprefix, int altname, int argc, char **argv, plotoptions_def plotoptions, verbose_definition verbose)
{
  char filename[MaxFilenameLength], txt[MaxFilenameLength];
  datafile_definition clone;
  if(altname) {
    strcpy(filename, argv[altname]);
  }else {
    if(extprefix == 0) {
      sprintf(txt, "%d.2dfs", file_number);
    }
    if(change_filename_extension(argv[argc-1], filename, txt, MaxFilenameLength, verbose) == 0) {
      return 0;
    }
  }
  if(verbose.verbose) {
    printf("Reading %s\n", filename);
  }
  closePSRData(&(twodfs_allinfo->datafile), 0, verbose);
  if(!openPSRData(&(twodfs_allinfo->datafile), filename, 0, 0, 1, 0, verbose)) {
    return 0;
  }
  if(PSRDataHeader_parse_commandline(&(twodfs_allinfo->datafile), argc, argv, verbose) == 0) {
    return 0;
  }
  if(verbose.verbose) {
    printf("%ldx%ld points read from 2dfs\n", twodfs_allinfo->datafile.NrBins, twodfs_allinfo->datafile.NrSubints);
  }
  if(twodfs_allinfo->datafile.NrPols > 1) {
    if(preprocess_polselect(twodfs_allinfo->datafile, &clone, 0, verbose) == 0) {
      printerror(verbose.debug, "ERROR pspecFig: Error selecting polarization channel 0.");
      return 0;
    }
    swap_orig_clone(&(twodfs_allinfo->datafile), &clone, verbose);
  }
  long i;
  double max;
  if(plotoptions.normalise_spectra) {
    for(i = 0; i < twodfs_allinfo->datafile.NrSubints * twodfs_allinfo->datafile.NrBins; i++) {
      if(i == 0 || fabs(twodfs_allinfo->datafile.data[i]) > max) {
 max = fabs(twodfs_allinfo->datafile.data[i]);
      }
    }
    for(i = 0; i < twodfs_allinfo->datafile.NrSubints * twodfs_allinfo->datafile.NrBins; i++) {
      twodfs_allinfo->datafile.data[i] /= 0.999*max;
    }
  }
  if(twodfs_allinfo->f2_min == 0.0 && twodfs_allinfo->f2_min == 0.0) {
    if(twodfs_allinfo->datafile.xrangeset) {
      twodfs_allinfo->f2_min = twodfs_allinfo->datafile.xrange[0];
      twodfs_allinfo->f2_max = twodfs_allinfo->datafile.xrange[1];
    }else {
      twodfs_allinfo->f2_min = -AverageProfile.NrBins/2.0;
      twodfs_allinfo->f2_max = -AverageProfile.NrBins/2.0 + AverageProfile.NrBins*(twodfs_allinfo->datafile.NrBins-1.0)/(float)twodfs_allinfo->datafile.NrBins;
    }
    float f2_resolution;
    f2_resolution = (twodfs_allinfo->f2_max - twodfs_allinfo->f2_min)/(float)(twodfs_allinfo->datafile.NrBins-1);
    twodfs_allinfo->f2_min -= 0.5*f2_resolution;
    twodfs_allinfo->f2_max += 0.5*f2_resolution;
    if(plotoptions.doFlip) {
      float swap;
      swap = twodfs_allinfo->f2_min;
      twodfs_allinfo->f2_min = -twodfs_allinfo->f2_max;
      twodfs_allinfo->f2_max = -swap;
    }
  }
  twodfs_allinfo->f3_min = 0;
  twodfs_allinfo->f3_max = 0.5;
  if(twodfs_allinfo->SelectP3Integrate) {
    if(twodfs_allinfo->P3IntegrateHigh > twodfs_allinfo->f3_max) {
      twodfs_allinfo->P3IntegrateHigh = twodfs_allinfo->f3_max;
    }
  }
  if(twodfs_allinfo->f3_min < 0) {
    twodfs_allinfo->f3_min = 0;
  }
  return 1;
}

ext: c

max_stars_repo_name: weltevrede/psrsalsa

max_issues_repo_name: weltevrede/psrsalsa

max_forks_repo_name: weltevrede/psrsalsa


text:
/*! \file deleter.h
    \brief gsl_interp_accelとgsl_splineのデリータを宣言・定義したヘッダファイル

    Copyright © 2015 @dc1394 All Rights Reserved.
    This software is released under the BSD 2-Clause License.
*/

#ifndef _DELETER_H_
#define _DELETER_H_

#pragma once

#include <gsl/gsl_spline.h>     // for gsl_interp_accel, gsl_interp_accel_free, gsl_spline, gsl_spline_free

namespace freefallsolveeom {
    //! A lambda expression.
    /*!
        gsl_interp_accelへのポインタを解放するラムダ式
        \param acc gsl_interp_accelへのポインタ
    */
    static auto const gsl_interp_accel_deleter = [](auto acc) {
        gsl_interp_accel_free(acc);
    };

    //! A lambda expression.
    /*!
        gsl_splineへのポインタを解放するラムダ式
        \param spline gsl_splineへのポインタ
    */
    static auto const gsl_spline_deleter = [](auto spline) {
        gsl_spline_free(spline);
    };
}

#endif  // _DELETER_H_

ext: h


text:
/*
 * BSD 3-Clause License
 *
 * Copyright (c) 2018, mtezych
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 *
 * * Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 *
 * * Neither the name of the copyright holder nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#ifndef CL_CONTEXT
#define CL_CONTEXT

#ifdef __APPLE__
	#include <OpenCL/cl.h>
#else
	#include <CL/cl.h>
#endif

#include <cl/Device.h>

#include <util/util.h>

#include <gsl/span>

#include <vector>
#include <cstddef>

namespace cl
{
	struct Platform;

	struct Context
	{
		cl_context clContext;

		//
		// @note: The NotifyCallback must be thread-safe.
		//
		using NotifyCallback = void (CL_CALLBACK *)
		(
			const char* errorInfo,
			const void* privateInfo, size_t cb,
			void* userData
		);

		Context
		(
			const Platform&         platform,
			gsl::span<const Device> devices
		);

		Context
		(
			const Platform& platform,
			cl_device_type  deviceType
		);

		~Context ();

		Context (Context&& context);
		Context (const Context& context) = delete;

		Context& operator = (Context&& context);
		Context& operator = (const Context& context) = delete;

		enum class Info : cl_context_info
		{
			NumDevices     = CL_CONTEXT_NUM_DEVICES,
			ReferenceCount = CL_CONTEXT_REFERENCE_COUNT,
		};

		template <Info info>
		auto GetInfo() const
		{
			auto infoSize = std::size_t { 0 };
			auto result = clGetContextInfo
			(
				clContext, util::enum_cast(info),
				0, nullptr, &infoSize
			);
			assert(result == CL_SUCCESS);

			auto infoBytes = std::vector<std::byte>
			{
				infoSize, std::byte { 0x00 }
			};
			result = clGetContextInfo
			(
				clContext, util::enum_cast(info),
				infoBytes.size(), infoBytes.data(), nullptr
			);
			assert(result == CL_SUCCESS);

			return InfoResult<info>::FromBytes(infoBytes);
		}

	private:

		template <Info info>
		struct InfoResult;
	};

	template <>
	struct Context::InfoResult<Context::Info::NumDevices>
	{
		static cl_uint
		FromBytes (const std::vector<std::byte>& infoBytes);
	};

	template <>
	struct Context::InfoResult<Context::Info::ReferenceCount>
	{
		static cl_uint
		FromBytes (const std::vector<std::byte>& infoBytes);
	};
}

#endif

max_stars_repo_name: mtezych/cpp

max_issues_repo_name: mtezych/cpp

max_forks_repo_name: mtezych/cpp


text:
/* min/bracketing.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* bracketing.c -- find an initial bracketing interval for a function to minimize */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_machine.h>

#include "min.h"

int
gsl_min_find_bracket(gsl_function *f,double *minimum,double * f_minimum,
		     double * x_lower, double * f_lower,
                     double * x_upper, double * f_upper,
		     size_t eval_max)
{
  /* The three following variables must be declared volatile to avoid storage
     in extended precision registers available on some architecture. The code
     relies on the ability to compare double values. As the values will be
     store in regular memory, the extended precision will then be lost and
     values that are different in extended precision might have equal
     representation in double precision. This behavior might break the
     algorithm.
   */
  volatile double f_left = *f_lower;
  volatile double f_right = *f_upper;
  volatile double f_center;
  double x_left = *x_lower;
  double x_right= *x_upper;
  double x_center;
  const double golden = 0.3819660;	/* golden = (3 - sqrt(5))/2 */
  size_t nb_eval = 0;

  if (f_right >= f_left)
    {
      x_center = (x_right - x_left) * golden + x_left;
      nb_eval++;
      SAFE_FUNC_CALL (f, x_center, &f_center);
    }
  else
    {
      x_center = x_right ;
      f_center = f_right ;
      x_right = (x_center - x_left) / golden + x_left;
      nb_eval++;
      SAFE_FUNC_CALL (f, x_right, &f_right);
    }

  do
    {
      if (f_center < f_left )
	{
	  if (f_center < f_right)
	    {
	      *x_lower = x_left;
	      *x_upper = x_right;
	      *minimum = x_center;
	      *f_lower = f_left;
	      *f_upper = f_right;
	      *f_minimum = f_center;
/*	      gsl_ieee_printf_double (&f_left);
	      printf(" ");
	      gsl_ieee_printf_double (&f_center);
	      printf("\n");*/
	      return GSL_SUCCESS;
	    }
	  else if (f_center > f_right)
	    {
	      x_left = x_center;
	      f_left = f_center;
	      x_center = x_right;
	      f_center = f_right;
	      x_right = (x_center - x_left) / golden + x_left;
	      nb_eval++;
	      SAFE_FUNC_CALL (f, x_right, &f_right);
	    }
	  else /* f_center == f_right */
	    {
	      x_right = x_center;
	      f_right = f_center;
	      x_center = (x_right - x_left) * golden + x_left;
	      nb_eval++;
	      SAFE_FUNC_CALL (f, x_center, &f_center);
	    }
	}
      else /* f_center >= f_left */
	{
	  x_right = x_center;
	  f_right = f_center;
	  x_center = (x_right - x_left) * golden + x_left;
	  nb_eval++;
	  SAFE_FUNC_CALL (f, x_center, &f_center);
	}
    }
  while (nb_eval < eval_max
	 && (x_right - x_left) > GSL_SQRT_DBL_EPSILON * ( (x_right + x_left) * 0.5 ) + GSL_SQRT_DBL_EPSILON);
  *x_lower = x_left;
  *x_upper = x_right;
  *minimum = x_center;
  *f_lower = f_left;
  *f_upper = f_right;
  *f_minimum = f_center;
  return GSL_FAILURE;
}


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#ifndef _MODEL1_H
#define _MODEL1_H

#include <stdbool.h>
#include <gsl/gsl_rng.h>

/* number of model parameters */
extern int num_params;

extern int which_level_update;
extern double *limits;     // limits from dnest
extern int thistask, totaltask;

extern DNestFptrSet *fptrset_thismodel;

/* functions */
void from_prior_thismodel(void *model);
void print_particle_thismodel(FILE *fp, const void *model);
double log_likelihoods_cal_thismodel(const void *model);
double perturb_thismodel(void *model);
void restart_action_model1(int iflag);
#endif

max_stars_repo_name: LiyrAstroph/CDNest


text:
/**
 *
 * @file testing_ssygst.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:18 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_smain.h"

#undef COMPLEX
#define REAL

static int check_transformation(int, int, int, float*, float*, int, float*, int, float);
static int check_factorization(int, float*, float*, int, int , float);

int testing_ssygst(int argc, char **argv)
{
    /* Check for number of arguments*/
    if (argc != 3) {
        USAGE("HEGST", "N LDA LDB",
              "   - N    : size of the matrices A and B\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - LDB  : leading dimension of the matrix B\n");
        return -1;
    }

    float eps = LAPACKE_slamch_work('e');
    int    N     = atoi(argv[0]);
    int    LDA   = atoi(argv[1]);
    int    LDB   = atoi(argv[2]);
    int    info_transformation, info_factorization;
    int    i, u;
    int    LDAxN = LDA*N;
    int    LDBxN = LDB*N;

    float *A1    = (float *)malloc(LDAxN*sizeof(float));
    float *A2    = (float *)malloc(LDAxN*sizeof(float));
    float *B1    = (float *)malloc(LDBxN*sizeof(float));
    float *B2    = (float *)malloc(LDBxN*sizeof(float));
    float *Ainit = (float *)malloc(LDAxN*sizeof(float));
    float *Binit = (float *)malloc(LDBxN*sizeof(float));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)||(!Ainit)||(!Binit)){
        printf("Out of Memory \n ");
        return -2;
    }

    /*----------------------------------------------------------
    *  TESTING SSYGST
    */

    /* Initialize A1 and A2 */
    PLASMA_splgsy(0., N, A1, LDA, 5198);
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, A1, LDA, Ainit, LDA);

    /* Initialize B1 and B2 */
    PLASMA_splgsy((float)N, N, B1, LDB, 4231);
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, B1, LDB, Binit, LDB);

    printf("\n");
    printf("------ TESTS FOR PLASMA SSYGST ROUTINE -------  \n");
    printf("        Size of the Matrix %d by %d\n", N, N);
    printf("\n");
    printf(" The matrices A and B are randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");

    /*----------------------------------------------------------
     *  TESTING SSYGST
     */

    for (i=0; i<3; i++) {
        for (u=0; u<2; u++) {
            memcpy(A2, Ainit, LDAxN*sizeof(float));
            memcpy(B2, Binit, LDBxN*sizeof(float));

            PLASMA_spotrf(uplo[u], N, B2, LDB);
            PLASMA_ssygst(itype[i], uplo[u], N, A2, LDA, B2, LDB);

            /* Check the Cholesky factorization and the transformation */
            info_factorization = check_factorization(N, B1, B2, LDB, uplo[u], eps);
            info_transformation = check_transformation(itype[i], uplo[u], N, A1, A2, LDA, B2, LDB, eps);

            if ( (info_transformation == 0) && (info_factorization == 0) ) {
                printf("***************************************************\n");
                printf(" ---- TESTING SSYGST (%s, %s) ....... PASSED !\n", itypestr[i], uplostr[u]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING SSYGST (%s, %s) ... FAILED !\n", itypestr[i], uplostr[u]);
                printf("************************************************\n");
            }
        }
    }

    free(A1);
    free(A2);
    free(B1);
    free(B2);
    free(Ainit);
    free(Binit);

    return 0;
}
/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
static int check_factorization(int N, float *A1, float *A2, int LDA, int uplo, float eps)
{
    float alpha = 1.0;
    float Anorm, Rnorm, result;
    int info_factorization;
    int i,j;

    float *Residual = (float *)malloc(N*N*sizeof(float));
    float *L1       = (float *)malloc(N*N*sizeof(float));
    float *L2       = (float *)malloc(N*N*sizeof(float));
    float *work                 = (float *)malloc(N*sizeof(float));

    memset((void*)L1, 0, N*N*sizeof(float));
    memset((void*)L2, 0, N*N*sizeof(float));

    LAPACKE_slacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, N, A2, LDA, L1, N);
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, N, A2, LDA, L2, N);
    if (uplo == PlasmaUpper)
       cblas_strmm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    else
       cblas_strmm(CblasColMajor, CblasRight, (CBLAS_UPLO)uplo, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);

    /* Compute the Residual || A - L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Residual, N,   work);
    Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1,       LDA, work);

    result = Rnorm / ( Anorm * N * eps );
    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n", result);

    if ( isnan(result) || isinf(result) || (result > 60.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*------------------------------------------------------------
 *  Check the Transformation to standard eigenvalue problem
 */

static int check_transformation(int itype, int uplo, int N, float *A1, float *A2, int LDA, float *B2, int LDB, float eps)
{
    float alpha = 1.0;
    float Anorm, Rnorm, result;
    int info_transformation;
    int i, j;
    char *str;

    float *Residual = (float *)malloc(N*N*sizeof(float));
    float *Aorig    = (float *)malloc(N*N*sizeof(float));
    float *work                 = (float *)malloc(N*sizeof(float));

    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'a',                N, N, A1, LDA, Residual, N);
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, N, A2, LDA, Aorig,    N);

    /* Rebuild the symmetry of A2 */
    if (uplo == PlasmaLower) {
        for (j = 0; j < N; j++)
            for (i = j+1; i < N; i++)
                Aorig[j+i*N] = (Aorig[i+j*N]);
    } else {
        for (i = 0; i < N; i++)
            for (j = i+1; j < N; j++)
                Aorig[j+i*N] = (Aorig[i+j*N]);
    }

    if (itype == 1) {
        if (uplo == PlasmaLower) {
            str = "L*A2*L'";
            cblas_strmm(CblasColMajor, CblasLeft,  CblasLower, CblasNoTrans,   CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
            cblas_strmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
        }
        else{
            str = "U'*A2*U";
            cblas_strmm(CblasColMajor, CblasLeft,  CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
            cblas_strmm(CblasColMajor, CblasRight, CblasUpper, CblasNoTrans,   CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
        }
    }
    else {
        if (uplo == PlasmaLower) {
            str = "inv(L')*A2*inv(L)";
            cblas_strsm(CblasColMajor, CblasLeft,  CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
            cblas_strsm(CblasColMajor, CblasRight, CblasLower, CblasNoTrans,   CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
        }
        else{
            str = "inv(U)*A2*inv(U')";
            cblas_strsm(CblasColMajor, CblasLeft,  CblasUpper, CblasNoTrans,   CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
            cblas_strsm(CblasColMajor, CblasRight, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);

        }
    }

    /* Compute the Residual ( A1 - W ) */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
            Residual[j*N+i] = Aorig[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Residual, N,   work);
    Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1,       LDA, work);

    result = Rnorm / (Anorm * N *eps);
    printf("============\n");
    printf("Checking the global transformation \n");
    printf("-- ||A1-%s||_oo/(||A1||_oo.N.eps) = %e \n", str, result );

    if (isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Transformation is suspicious ! \n");
        info_transformation = 1;
    }
    else {
        printf("-- Transformation is CORRECT ! \n");
        info_transformation = 0;
    }

    free(Residual); free(Aorig); free(work);

    return info_transformation;
}