#include "../cpgplot.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define C 299792458
#define DAY 86400
#define PI 3.141592653589793
#define E 2.718281828
#define GMSUN 1.32712438E+20
#define TSUN 0.000004925490947
#define RSUN 6.96E+8
#define LIM 1E-10
#define DE 1E-5
#define DM1 0.01
#define NM1 1000 
#define NHH 80
#define NHV 50
#define NLMAX 10
#define NM 1024 
#define NM2 1048576 
#define NCONT 3
#define one3rd 0.33333333333
#define two3rd 0.66666666666
#define four3rd 1.33333333333
#define five3rd 1.66666666666
#define dsin 0.0001

 /* ****************************************************************************************************** */

float fm(float x, float s, float f, float m1)
{
  float r;
  r = x*x*x*s*s*s - f*x*x - 2*m1*f*x - f*m1*m1;
  return r;
}
 
float gm(float x, float gamma, float gammak, float m1)
{
  float r;
  r = gamma * pow((m1 + x), four3rd) - gammak*x*(m1 + 2*x);
  return r;
}
 
float mass(int k, float mi, float a1, float a2, float a3)
{
  int i;
  float m, xn, xo, dif;
 
  /* This routine finds zero of previously declared functions */

  xo = 0.0;
  xn = mi;
  i = 0;
  dif = sqrt((xn - xo)*(xn-xo));
  while (dif > LIM)
    {
      if (k ==0) m = ( fm(xn + 0.5*DE, a1, a2, a3) - fm(xn - 0.5*DE, a1, a2, a3) )/DE;
      if (k ==1) m = ( gm(xn + 0.5*DE, a1, a2, a3) - gm(xn - 0.5*DE, a1, a2, a3) )/DE;
      xo = xn;
      if (k ==0) xn = xo - fm(xo, a1, a2, a3)/m;
      if (k ==1) xn = xo - gm(xo, a1, a2, a3)/m;

      dif = sqrt((xn - xo)*(xn-xo));
      i++;
      if (i > 100)
        dif = 0;
    }
  return xn;
}

main( )
{

  /* ****************************************************************************************************** */

  /* PRELIMINARY STAGE: VARIABLE DECLARATIONS */

  /* declare orbital and derived parameters */
  double f, ftest, m1, m2, inc, s, ds, asini, pb, n, e;
  double wdot, ewdot, odot, dodot, M, Ml, Mu, Mp;
  double gamma, gammam, egamma, gammak, pbdot, pbdotm, epbdot, pbdotk;
  double rs, rsm, ers, ss, ssm, ess, R, Rm, eR;
  
  float q, a, rl, rc, mue2, w, exponent;
  double mtotal, mptotal, mctotal, total, height, height0;
  double cosi, sini, sin1, sin2, cos1, cos2, cosidiff, sinidiff;

  /* declare files for p.d.f.s */
  FILE *PM;
  FILE *CM;
  FILE *FGAMMA;
  FILE *FPBDOT;
  FILE *FRS;
  FILE *FSS;
  FILE *FR;

  /* declare counters and flags*/
  int b, c, d1, d2, d3, i, j, k, l;
  int s0, s1, s2, s3, s4, s5, t;
  int pk0, pk1, pk2, pk3, pk4, pk5, pk6;
  int pnow, pbefore, finish;

  /* declare graphical parameters */
  float dm, dm1, dm2, m1l, m1u, m2l, m2u, c1, c2, mlow, mhigh;

  /* declare line and probability/countour parameters */
  int ni, nilines, nm, nmlines, ng, nglines, npb, npblines;
  int  nrs, nrslines, nss, nsslines, nR, nRlines;
  float si[NLMAX], sigmam[NLMAX], sigmag[NLMAX], sigmapb[NLMAX];
  float sigmars[NLMAX], sigmass[NLMAX], sigmaR[NLMAX];
  float tr[6], alev[NCONT], alevc[NCONT];
  float mpalev[NCONT], mcalev[NCONT], mpalevy[NCONT], mcalevy[NCONT];
  float target[NCONT], mtarget[NCONT];
  float mx[NM], my[NM], nx[2], ny[2];
  float fnp[NM], fnc[NM], fmp[NM], fmc[NM], np[NM], nc[NM];
  float fnpc[NM2], npc[NM2], npcmax;
  double fgamma[NM], fpbdot[NM], frs[NM], fss[NM], fR[NM];
  double npmax, ncmax, flagp[NCONT], flagc[NCONT];

  /* ****************************************************************************************************** */

  /* Very basic level: calculating the mass function if user does not know it */

  s0 = 'n';  s1 = 'n'; s2 = 'n', s3 = 'n', s4 = 'n', s5 = 'n', t = 0;

  printf("\n For all questions below, the t option is to terminate the program without further questions.\n");
  printf("\n Do you know the mass function? (y/n/t)\n");
  s0 = getchar();

  if (s0 == 'y')
    {
      printf("\n Input mass function (in solar masses): \n");
      scanf("%lf", &f);
      s1 = getchar();
    }
  if (s0 == 'n')
    {
      printf("\n Input x (light seconds): \n");
      scanf("%lf", &asini);
      printf("\n Input the orbital period (days): \n");
      scanf("%lf", &pb);
      n = 2 * PI / (pb * DAY);
      /* Now: Calculate f in solar masses... */
      f = pow((asini * C), 3) * pow(( 2 * PI /(pb * DAY)), 2) / GMSUN;
      printf ("\n The mass function is %e solar masses.\n", f);
      c = getchar();
    }
  if (s0 == 't')
    t = 1;

  /* ****************************************************************************************************** */

  /* First level of real use of the program: companion mass calculator */

  if (t == 0)
    {
      printf("\n Do you want to calculate the companion mass assuming a pulsar mass and sin i? (y/n/t)\n");
      s1 = getchar();
      
      if (s1 == 'y')
	{
	  printf("\n What is the assumed mass of neutron star: \n"); 
	  scanf("%lf", &m1);
	  printf("\n What is the inclination? (degrees): \n");
	  scanf("%lf", &inc);
	  s = sin(inc * PI/180.0);
	  
	  m2 = m1*m1*f/(s*s*s); m2 = pow(m2,(1/3));
	  m2 = (m1+m2)*(m1+m2)*f/(s*s*s); m2 = pow(m2,(1/3));
	  
	  m2 = mass(0,m2,s,f,m1);
	  
	  printf("\n Mass of companion: %e solar masses.\n", m2);
	  
	  printf("\n Do you want to calculate geometric parameters? (y/n/t)\n");
	  
	  s2 = getchar();
	  s2 = getchar();
	  
	  if ( s2 == 'y' )
	    {
	  /* OK, let's calculate the geometric parameters */
	      
	  /* Do we know all we need? If we only know the mass function, we need to find out asini, PB*/
	      
	      if (s0 == 'y')
		{
		  printf("\n Input the orbital period (days): \n");
		  scanf("%lf", &pb);
		  n = 2 * PI / (pb * DAY);
		  asini = f * TSUN / (n*n);
		  asini = pow (asini, one3rd);
		}
	      else
		{
		  printf("\n We will use the previously introduced orbital period, %f days.", pb);
		}
	      
	      /* We will also need the orbital eccentricity at some point here */
	      
	      printf("\n Input the orbital eccentricity: \n", q);
	      scanf("%lf", &e);
	      
	      /* Let's now print the other parameters:
		 - mass quocient
		 - separation
		 - Roche Lobe Radius
		 - radius of companion object if it is a white dwarf */ 
	      q = m2 / m1;
	      printf("\n Mass of companion / mass of pulsar: %e.\n", q);
	      a = asini * C * (1 + m1/m2) / s;
	      printf("\n Mean separation: %e m (%f solar radii)\n", a, a/RSUN);
	      printf("\n Separation at periastron: %e m (%f solar radii)\n", a*(1-e), a*(1-e)/RSUN);
	      printf("\n Separation at apastron: %e m (%f solar radii)\n", a*(1+e), a*(1+e)/RSUN);
	      rl = a* 0.49 * pow (q, two3rd) / (0.6 * pow(q, two3rd) + log(1 + pow(q, one3rd))/log(E));
	      printf("\n Roche lobe radius for circular orbit: %e m (%f solar radii)\n", rl, rl/RSUN);
	      mue2 = 1;
	      rc = m2 * pow(mue2,5) / 0.7011;
	      rc = 1E+7*pow(rc,-one3rd);
	      printf("\n Companion radius, m = 2: %e m (%f solar radii)\n", rc, rc/RSUN); 
	      printf("\n Fraction of power from the pulsar intercepted by the companion: %f\n", rc*rc/(4*a*a));
	      w = m2 * GMSUN / rc;
	      printf("\n Gravitational potential at the surface is %e (m/s)^2\n", w);
	      mue2 = 0.5;
	      rc = m2 * pow(mue2,5) / 0.7011;
	      rc = 1E+7*pow(rc,-one3rd);
	      printf("\n Companion radius, m = 1: %e m (%f solar radii)\n", rc, rc/RSUN);
	      printf("\n Fraction of power from the pulsar intercepted by the companion: %f\n", rc*rc/(4*a*a));
	      w = m2 * GMSUN / rc;
	      printf("\n Gravitational potential at the surface is %e (m/s)^2\n", w);
	    }
	  if (s2 == 't')
	    t = 1;
	  
	  if (t == 0)
	    {
	      printf("\n Do you want a prediction of the relativistic effects for these masses and inclination? (y/n/t) \n");
	      s3 = getchar();
	      s3 = getchar();
	      
	      if (s3 == 'y')
		{
		  
		  /* Do we know everything we need? */
		  /* If the answer to s2 is yes, we know everything, inclusing the eccentricity. */
		  /* If not, then we have at least to find out what the eccentricity is. */
		  
		  if (s2 == 'n')
		    {
		      printf("\n Input the orbital eccentricity: \n");
		      scanf("%lf", &e);
		      
		      /* If the answer to the first question was 'yes', then we are still missing the orbital period */
		      
		      if (s0 == 'y')
			{
			  printf("\n Input the orbital period (days): \n");
			  scanf("%lf", &pb);
			  n = 2 * PI / (pb * DAY);
			}
		      else
			{
			  printf("\n We will use the previously introduced orbital period, %f days", pb);
			}
		    }
		  else
		    {
		      printf("\n We will use the previously introduced orbital eccentricity, %f", e);
		    }
		  
		  /* now we should know everything we need... */
		  M = m1 + m2;
		  
		  /* calculate the expected omegadot */
		  wdot = pow ( n , five3rd) / (1 - e * e) * pow( (M*TSUN) , two3rd) * 3;
		  wdot = wdot * 180 * 365.2425 * DAY /PI;
		  printf("\n The periastron advance is %f degrees/yr.\n", wdot);
		  gammak = e * pow( n , -one3rd )* pow( TSUN , two3rd );
		  gamma = gammak * m2 * ( M + m2 ) * pow ( M, -four3rd );
		  printf("\n The Einstein delay (gamma) is %f seconds.\n", gamma);
		  pbdotk = -192 * PI * pow ( (n * TSUN), five3rd) / 5.0;
		  pbdotk = pbdotk * (1 + e*e*73.0/24.0 + 37.0*pow(e , 4)/96.0);
		  pbdotk = pbdotk * pow((1 - e*e),-3.5);
		  pbdot = pbdotk * m1 * m2 / pow(M, one3rd);
		  printf("\n The orbital decay is %e.\n", pbdot);
		  c = getchar();
		}
	      if (s3 == 't')
		t = 1;
	    }
	}
      if (s1 == 't')
	t = 1;
      else
	s4 = getchar();
    }

  /* ****************************************************************************************************** */

  /* Higher level of use: estimate masses if one has an omega-dot */

  if (t == 0)
    {
      printf("\n Do you want to estimate masses given a periastron advance? (y/n/t)\n ");
      s4 = getchar();

      if (s4 == 'y')
	{

	  /* OK: same conversation: do we know everything we need? */
	  /* Here we have to contend with the fact that the answer to s1 might have been negative */
	  
	  /* If the answer to s2 is yes, we know everything, inclusing the eccentricity. */
	  /* If not, then we have at least to find out what the eccentricity is. */
	  
	  
	  /* If the answer to s3 was 'yes', we know e for sure */

	  if ((s3 == 'n') && (s2 == 'n'))
	    {
	      printf("\n Input the orbital eccentricity: \n");
	      scanf("%lf", &e);
	      
	      /* Only if the answer to s2 e s3 was 'no' and the answer to the first question was 'yes' can we be still missing the orbital period */
	      
	      if (s0 == 'y')
		{
		  printf("\n Input the orbital period (days): \n");
		  scanf("%lf", &pb);
		  n = 2 * PI / (pb * DAY);
		}
	      else
		{
		  printf("\n We will use the previously introduced orbital period, %f days.", pb);
		}
	      /* The nice thing here is that we don't need the mass and sin i assumptions. So, with "n" and "e", we have everything we need */
	    }
	  else
	    {
	      printf("\n We will use the previously introduced orbital eccentricity, %f and oprbital period, %f days", e, pb);
	    }
	  
	  printf("\n Input the periastron advance (degrees/year): \n");
	  scanf("%lf", &wdot);
	  
	  printf("\n What is the uncertainty in the periastron advance?\n(one sigma, degrees/year)\n");
	  scanf("%lf", &ewdot);
	  
	  /* NOW: DO CALCULATIONS */
	  
	  wdot = wdot * PI / (180 * 365.2425 * DAY);
      
	  ewdot = ewdot * PI / (180 * 365.2425 * DAY);
	  
	  M = wdot * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  M = pow (M,1.5);
	  printf("\n The total mass is %f solar masses.\n", M);

	  Ml = (wdot - 1 * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  Ml = pow (Ml,1.5);
	  Mu = (wdot + 1 * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  Mu = pow (Mu,1.5);
	  
	  printf("\n Lower limit (1-sigma) is %f solar masses.\n", Ml);
	  printf("\n Upper limit (1-sigma) is %f solar masses.\n", Mu);

	  Ml = (wdot - 2 * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  Ml = pow (Ml,1.5);
	  Mu = (wdot + 2 * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  Mu = pow (Mu,1.5);

	  printf("\n Lower limit (2-sigma) is %f solar masses.\n", Ml);
	  printf("\n Upper limit (2-sigma) is %f solar masses.\n", Mu);

	  Ml = (wdot - 3 * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  Ml = pow (Ml,1.5);
	  Mu = (wdot + 3 * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
	  Mu = pow (Mu,1.5);
	  
	  printf("\n Lower limit (3-sigma) is %f solar masses.\n", Ml);
	  printf("\n Upper limit (3-sigma) is %f solar masses.\n", Mu);
	  
	  c = 'y';
	  s = 1;
	  
	  while (c == 'y')
	    {
	      
	      m2 = M * M * f * pow (s , -3);
	      m2 = pow( m2, one3rd);
	      
	      printf("\n For sin i = %lf the mass of the companion is %f and the mass of the pulsar is %f.\n", s, m2, M-m2);
	      
	      m2 = Ml * Ml * f * pow (s , -3);
	      m2 = pow( m2, one3rd);
	      
	      printf("\n In the 3-sigma lower limit of total mass, the mass of the companion is %f .\n", m2);
	      printf("   the mass of the pulsar is %f .\n", Ml-m2);
	      
	      m2 = Mu * Mu * f * pow (s , -3);
	      m2 = pow( m2, one3rd);
	      
	      printf("\n In the 3-sigma upper limit of total mass, the mass of the companion is %f .\n", m2);
	      printf("   the mass of the pulsar is %f .\n", Mu-m2);
	      
	      printf ("\n Do you want to calculate the pulsar and companion masses for another inclination? (y/n)\n");
	      
	      c = getchar();
	      c = getchar();
	      
	      if (c == 'y')
		{
		  printf("\n What is the inclination? (degrees): \n");
		  scanf("%lf", &inc);
		  s = sin(inc * PI/180.0);
		}
	    }
	}
      if (s4 == 't')
	t = 1;
    }

  /* New version comment:
     Omega-dot cycle ended here! Knowing omega-dot is no longer a necessary condition to make the plot! */


  /* ****************************************************************************************************** */

  /* Even higher level of use: make graphics */

  if (t == 0)
    {
      printf("\n Do you want to make a mass-mass diagram? (y/n)\n");
      s5 = getchar();
      s5 = getchar();
      
      if (s5 == 'y')
	{
	  
	  /* Question number one: do we know the mass function? Yes - does not matter what answer to s0 was. */
	  
	  /* Because we know the mass function, we can draw the constant inclination lines. Pb and e don't matter. */
	  
	  printf("\n How many constant inclination lines do you want to print?\n");
	  scanf("%d", &nilines);
	  si[0] = 1;
	  for (ni = 1; ni < nilines; ni++)
	    {
	      printf("\n What is the inclination of line number %d?\n", ni+1);
	      scanf("%lf", &inc);
	      si[ni] = sin(inc * PI/180.0);
	      printf("\n The sine of %f degrees is %f\n", inc, si[ni]);
	    }
	  
	  /* Question number two: do we know the precession of periastron? */
	  /* If s4 = 'y', we're fine, and we surely know PB and e and om-dot. Otherwise, we have to ask */
	  
	  if (s4 == 'n')
	    {
	      if ((s3 == 'n') && (s2 == 'n'))
		{
		  printf("\n Input the orbital eccentricity: \n");
		  scanf("%lf", &e);
		  
		  /* Only if the answer to s2 e s3 was 'no' and the answer to the first question was 'yes' can we be still missing the orbital period */
	      
		  if (s0 == 'y')
		    {
		      printf("\n Input the orbital period (days): \n");
		      scanf("%lf", &pb);
		      n = 2 * PI / (pb * DAY);
		    }
		  else
		    {
		      printf("\n We will use the previously introduced orbital period, %f days,", pb);
		    }
		}
	      else
		{
		  printf("\n We will use the previously introduced orbital eccentricity, %f.", e);
		}
	      
	      printf("\n Do you know the periastron advance? (y/n)\n");
	      pk1 = getchar();
	      pk1 = getchar();
	      
	      if (pk1 == 'y')
		{
		  printf("\n Input the periastron advance (degrees/year): \n");
		  scanf("%lf", &wdot);
		  
		  printf("\n What is the uncertainty in the periastron advance?\n(one sigma, degrees/year)\n");
		  scanf("%lf", &ewdot);
		  
		  wdot = wdot * PI / (180 * 365.2425 * DAY);
		  ewdot = ewdot * PI / (180 * 365.2425 * DAY);
		}
	    }
	  
	  else
	    {
	      pk1 = 'y';
	      printf("\n We will use the previously introduced orbital period, eccentricity and omega-dot");
	    }
	  
	  if (pk1 == 'y')
	    {
	      printf("\n How many constant total mass lines do you want to print?\n");
	      scanf("%d", &nmlines);
	      sigmam[0] = 0;
	      for (nm = 1; nm < nmlines; nm++)
		{
		  printf("\n How many sigmas is line %d above measurement?\n", nm+1);
		  scanf("%f", &sigmam[nm]);
		}
	    };
	  
	  printf("\n Do you know the Einstein delay (gamma)? (y/n)\n ");
	  pk2 = getchar();
	  pk2 = getchar();

	  if (pk2 == 'y') 
	    {
	      printf("\n Input gamma (seconds): \n");
	      scanf("%lf", &gammam);
	      printf("\n What is the uncertainty of gamma (one sigma, in seconds)?\n");
	      scanf("%lf", &egamma);      
	      gammak = e*pow(n,-one3rd)*pow(TSUN,two3rd);
	      
	      printf("\n How many constant gamma lines do you want to print?\n");
	      scanf("%d", &nglines);
	  
	      sigmag[0] = 0;
	      for (ng = 1; ng < nglines; ng++)
		{
		  printf("\n How many sigmas is line %d above measurement?\n", ng+1);
		  scanf("%f", &sigmag[ng]);
		}
	    };
	  
	  printf("\n Do you know the orbital decay (PBdot)? (y/n)\n ");
	  pk3 = getchar();
	  pk3 = getchar();
	  
	  if (pk3 == 'y') 
	    {
	      printf("\n Input PBdot (divided by 10^-12): \n");
	      scanf("%lf", &pbdotm); pbdotm = pbdotm * 1E-12;
	      printf("\n What is the uncertainty of PBdot (one sigma, in 10^-12)?\n");
	      scanf("%lf", &epbdot); epbdot = epbdot * 1E-12;
	      
	      printf("\n How many constant pbdot lines do you want to print?\n");
	      
	      scanf("%d", &npblines);
	      
	      sigmapb[0] = 0;
	      for (npb = 1; npb < npblines; npb++)
		{
		  printf("\n How many sigmas is line %d above measurement?\n", npb+1);
		  scanf("%f", &sigmapb[npb]);
		};
	      
	      /* Let us calculate pbdotk, as above, since we _must_ know n and e by now */
	      /* (this might be a repeat calculation, but I don't care - it does not take that long */
	      
	      pbdotk = -192 * PI * pow ((n * TSUN), five3rd) / 5.0;
	      pbdotk = pbdotk * (1 + e*e*73.0/24.0 +  pow(e , 4)*37.0/96.0);
	      pbdotk = pbdotk * pow((1 - e*e),-3.5);
	    };
	  
	  printf("\n Do you know the Shapiro parameter r? (same as m2 if GR is true, y/n)\n ");
	  pk4 = getchar();
	  pk4 = getchar();
	  
	  if (pk4 == 'y') 
	    {
	      printf("\n Input r (in solar masses): \n");
	      scanf("%lf", &rsm);
	      printf("\n What is the uncertainty of this parameter?\n");
	      scanf("%lf", &ers);
	      
	      printf("\n How many constant r lines do you want to print?\n");
	      
	      scanf("%d", &nrslines);
	      
	      sigmars[0] = 0;
	      for (nrs = 1; nrs < nrslines; nrs++)
		{
		  printf("\n How many sigmas is line %d above measurement?\n", nrs+1);
		  scanf("%f", &sigmars[nrs]);
		}
	    }
	  
	  printf("\n Do you know the Shapiro parameter s (same as sin i if GR is true)? (y/n)\n ");
	  pk5 = getchar();
	  pk5 = getchar();
      
	  if (pk5 == 'y') 
	    {
	      printf("\n Input s: \n");
	      scanf("%lf", &ssm);
	      printf("\n What is the uncertainty of this parameter?\n");
	      scanf("%lf", &ess);
	      
	      printf("\n How many constant s lines do you want to print?\n");
	      
	      scanf("%d", &nsslines);
	      
	      sigmass[0] = 0;
	      for (nss = 1; nss < nsslines; nss++)
		{
		  printf("\n How many sigmas is line %d above measurement?\n", nss+1);
		  scanf("%f", &sigmass[nss]);
		}
	    }

	  printf("\n Do you know the mass ratio? (y/n)\n ");
	  pk6 = getchar();
	  pk6 = getchar();

	  if (pk6 == 'y') 
	    {
	      printf("\n Input R: \n");
	      scanf("%lf", &Rm);
	      printf("\n What is the uncertainty of this parameter?\n");
	      scanf("%lf", &eR);
	      
	      printf("\n How many constant R lines do you want to print?\n");
	      
	      scanf("%d", &nRlines);
	      
	      sigmaR[0] = 0;
	      for (nR = 1; nR < nRlines; nR++)
		{
		  printf("\n How many sigmas is line %d above measurement?\n", nR+1);
		  scanf("%f", &sigmaR[nR]);
		}
	    };
	  
	  printf("\n Please input now the parameters related to the plot.\n");
	  
	  
	  printf("\n Figure details:\n");
	  printf("\n Do you want the region below the mass function to be hatched? (y/n)\n ");
	  d1 = getchar();
	  d1 = getchar();
	  
	  printf("\n Do you want to have a contour plot? (y/n)\n ");
	  d2 = getchar();
	  d2 = getchar();
	  
	  printf("\n Do you want to have the windows for the 1-d p.d.f.s? (y/n)\n ");
	  d3 = getchar();
	  d3 = getchar();
	  
	  /* Ask about output, make display, etc */
	  
	  printf("\n Please indicate lower and upper limits for X axis (in solar masses). \n"); 
	  scanf("%f", &m1l);
	  scanf("%f", &m1u);
	  
	  printf("\n Please indicate lower and upper limits for Y axis (in solar masses). \n"); 
	  scanf("%f", &m2l);
	  scanf("%f", &m2u);
	  
	  c = getchar();
	  
	  /* ************************************************************************* */
	  
	  /* NOW, THE FIGURE: START WITH BASIC LINES AND HATCHING */
	  
	  if(cpgbeg(0, "?", 1, 1) != 1)
	    return EXIT_FAILURE;
	  
	  cpgsch(1.4);
	  if (d3 == 'y')
	    {
	      cpgvsiz(1.0, 5.5, 1.0, 5.5);
	    }
	  else
	    {
	      cpgvsiz(1.0, 7.5, 1.0, 7.5);
	    }
	  
	  cpgswin(m1l, m1u, m2l, m2u);
	  cpglab("Pulsar mass (M\\d\\(2281)\\u)", "Companion mass (M\\d\\(2281)\\u)", "");
	  
	  dm = (m1u - m1l)/(NM1-1);	  
	  
	  /* print band with previously measured masses */
	  
	  cpgsci(15);
	  cpgsfs(1);
	  
	  mlow = 1.20;
	  mhigh = 1.4408;
	  
	  mx[0] = mlow;
	  my[0] = m2l;
	  mx[1] = mlow;
	  my[1] = m2u;
	  mx[3] = mhigh;
	  my[3] = m2l;
	  mx[2] = mhigh;
	  my[2] = m2u;

	  cpgpoly(4, mx, my);
	  
	  mx[0] = m1l;
	  my[0] = mlow;
	  mx[1] = m1l;
	  my[1] = mhigh;
	  mx[2] = m1u;
	  my[2] = mhigh;
	  mx[3] = m1u;
	  my[3] = mlow;
	  
	  cpgpoly(4, mx, my);
	  
	  cpgsci(1);
	  
	  /* print m2 assuming sin i = x and the known mass function, f */  
	  
	  cpgslw(5);
	  /* First round: print sin i = 1, then do other rounds */
	  
	  /* I found a really cool alternative way of calculating the constant inclination curves */
	  
	  for (ni = 0; ni < nilines; ni++)
	    {
	      s = si[ni];
	      /* calculate m2 for mp = 0 */
	      m2 = f/(s*s*s);
	      mx[0] = 0; my[0] = m2;
	      
	      for (k = 1; k < NM1; k++)
		{
		  /* First: cycle through mass ratios */
		  R = 10.0/(1.0*k);
		  mx[k] = m2 * (1+R)*(1+R)/(R*R*R);
		  my[k] = R*mx[k];
		}
	      cpgline(NM1, mx, my);
	      cpgslw(2);  /* following lines to be drawn with this width */
	    };

	  /* Second round: print total mass lines - only if pk1 = 'y'*/
	  
	  cpgsls(2); /*remaining lines, if any, are done in style 2 */
	  cpgslw(1);
	  
	  /* The measured omega line is done automatically */
	  
	  if (pk1 == 'y')
	    {
	      for (nm = 0; nm < nmlines; nm++)
		{
		  Mp = (wdot + sigmam[nm] * ewdot) * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
		  Mp = pow (Mp,1.5);

		  /* How does the total mass translate at the mass function (sin i = 1)? */
		  my[0] = pow( ( f * Mp * Mp ), one3rd);
		  mx[0] = Mp - my[0];
	      
		  mx[1] = 0;
		  my[1] = Mp;
		  
		  cpgline(2, mx, my);
		}
	    };
	  
	  cpgsls(3); /*remaining lines, if any, are done in style 3 */
	  cpgslw(1);
	  
	  /* Now, the gamma lines - only if pk2 = 'y' */
	  
	  if (pk2 == 'y')
	    {
	      for (ng = 0; ng < nglines; ng++)
		{
		  gamma = gammam + sigmag[ng]*egamma;
		  
		  /* initialize the value of m2 for m1 = 0 */
		  m2=gamma/(2*gammak);
		  m2=pow(m2,1.5);
		  
		  /* printf("\n e = %f, n = %f\n", e, n);
		     printf("\n Initial mass guess = %f solar masses\n", m2);*/
		  if (m2 < m2u)
		    {
		      for (k = 0; k < NM1; k++)
			{
			  m1 = k * dm + m1l;
			  mx[k] = m1;
			  m2 = mass(1,m2,gamma,gammak,m1);
			  my[k] = m2;
			  /* see if, for sin i = 1, we get a value larger than the mass function */
			}  
		      cpgline(NM1, mx, my);
		    }
		}
	    }
	  
	  cpgsls(2); /*remaining lines, if any, are done in style 3 */
	  cpgslw(1);
	  
	  /* now, the pbdot lines - only if pk3 = 'y' */
	  
	  if ( pk3 == 'y' )
	    {
	      for (npb = 0; npb < npblines; npb++)
		{
		  pbdot = pbdotm + sigmapb[npb]*epbdot;
		  
		  /* ludicrous values for first entry in vector - not likely to be displayed */
		  mx[0] = 100;
		  my[0] = 0;
		  for (k = 1; k < NM1; k++)
		    {
		      /* First: cycle through mass rations */
		      R = 0.1 * k;
		      m2 = pbdot * pow((1+R), one3rd)/(R * pbdotk);
		      m2 = pow(m2,0.6);
		      m1 = R * m2;
		      mx[k] = m1;
		      my[k] = m2;
		    }
		  cpgline(NM1, mx, my);
		}
	    }
	  
	  /*print constant mass lines - easy :-) */
	  
	  cpgsls(4); /*remaining lines, if any, are done in style 4 */
	  
	  if ( pk4 == 'y' )
	    {
	      for (nrs = 0; nrs < nrslines; nrs++)
		{
		  rs = rsm + sigmars[nrs]*ers;

		  mx[0] = m1l;
		  my[0] = rs;
		  
		  mx[1] = m1u;
		  my[1] = rs;
		  
		  cpgline(2, mx, my);
		}
	    }
	  
	  /* print s lines - sames as constant inclination lines */
	  
	  if ( pk5 == 'y' )
	    {
	      for (nss = 0; nss < nsslines; nss++)
		{
		  ss = ssm + sigmass[nss]*ess;
		  
		  /* calculate m2 for mp = 0 */
		  m2 = f/(ss*ss*ss);
		  mx[0] = 0; my[0] = m2;
		  
		  for (k = 1; k < NM1; k++)
		    {
		      /* First: cycle through mass rations */
		      R = 10.0/(1.0*k);
		      mx[k] = m2 * (1+R)*(1+R)/(R*R*R);
		      my[k] = R*mx[k];
		    }
		  cpgline(NM1, mx, my);
		}
	    }
	  
	  if ( pk6 == 'y' )
	    {
	      for (nR = 0; nR < nRlines; nR++)
		{
		  R = Rm + sigmaR[nR]*eR;
		  
		  mx[0] = 0;
		  my[0] = 0;
		  
		  mx[1] = m1u;
		  my[1] = m1u/R;
		  
		  cpgline(2, mx, my);
		}
	    }
	  
	  /* now : put the box around the plot */
	  
	  cpgsls(1);
	  cpgslw(1);
	  cpgbox("BCTNS", 0.0, 0, "BCTNS", 0.0, 0);
	  
	  /* Do the hatching? */
	  
	  if (d1 == 'y')
	    {
	      cpgslw(1);
	      cpgsls(1);
	      ds = (m1u - m1l) / NHH;
	      
	      for (k = 0; k < NHH; k++)
		{

		  /* now, we specify points near the X axis */
		  mx[0] = m1l + k * ds;
		  my[0] = m2l;

		  M = mx[0] + my[0];
		  /* How does the total mass translate at the mass function (sin i = 1)? */
		  my[1] = pow( ( f * M * M ), one3rd);
		  mx[1] = M - my[1];
	      
		  cpgline(2, mx, my);
		}      
	      
	      for (k = 0; k < NHH; k++)
		{
		  
		  /* now, we specify points near the Y axis */
		  mx[0] = m1u;
		  my[0] = m2l + k * ds;
		  
		  M = mx[0] + my[0];
		  /* How does the total mass translate at the mass function (sin i = 1)? */
		  my[1] = pow( ( f * M * M ), one3rd);
		  mx[1] = M - my[1];
	      
		  cpgline(2, mx, my);
		}
	    }
	  
	  /* ************************************************************************* */
	  
	  /* Do the contours? */
	  /* Calculate probability density functions? */
	  
	  /* In any of these eventualities, we have to calculate 2-d p.d.f. */
	  
	  if ((d2 == 'y')||(d3 == 'y'))
	    {
	      
	      /* go through a cycle - ask what quantities will be used to calculate the pdf */
	      
	      
	      printf("\n Will the inclinations be used in calculating the pdf? (y/n)\n");
	      
	      pk0 = getchar();
	      
	      if (pk1 == 'y')
		{
		  printf("\n Will the omega-dot be used in calculating the pdf? (y/n)\n");
		  
		  pk1 = getchar();
		  pk1 = getchar();
		}	  	
	      
	      
	      if (pk2 == 'y')
		{
		  printf("\n Will gamma be used in calculating the pdf? (y/n) or do you want a prediction (p)?\n");
		  
		  pk2 = getchar();
		  pk2 = getchar();
		}	  	
	      

	      if (pk3 == 'y')
		{
		  printf("\n Will the pb-dot be used in calculating the pdf (y/n) or do you want a prediction (p)?\n");
		  
		  pk3 = getchar();
		  pk3 = getchar();
		}	  	


	      if (pk4 == 'y')
		{
		  printf("\n Will the Shapiro r be used in calculating the pdf (y/n) or do you want a prediction (p)?\n");
		  
		  pk4 = getchar();
		  pk4 = getchar();
		}	  	
	      
	      
	      if (pk5 == 'y')
		{
		  printf("\n Will the Shapiro s be used in calculating the pdf (y/n) or do you want a prediction (p)?\n");
		  
		  pk5 = getchar();
		  pk5 = getchar();
		}	  	
	      
	      
	      if (pk6 == 'y')
		{
		  printf("\n Will the mass ratio be used in calculating the pdf (y/n) or do you want a prediction (p)?\n");
		  
		  pk6 = getchar();
		  pk6 = getchar();
		}	  	
	      
	      /* calculate conversion factors between mass and array number
		 these numbers will be needed in any case */
	      
	      c1 = (m1u - m1l)/(NM*1.0);
	      c2 = (m2u - m2l)/(NM*1.0);
	      
	      /* Initialize probability vectors */
	      
	      for (i = 0; i< NM; i++)
		{
		  np[i] = 0.0;
		  nc[i] = 0.0;
		}
	      for (i = 0; i< NM2; i++)
		npc[i] = 0.0;
	      
	      mtotal = 0;
	      
	      
	      /* If we have an omega-dot, then we assume that it is the most constraining factor - we go into hyper-precision mode */
	      
	      printf("\n Calculating pdf with ");
	      if ( pk0 == 'y' ) printf("inclinations, ");
	      if ( pk1 == 'y' ) printf("omega-dot, ");
	      if ( pk2 == 'y' ) printf("gamma, ");
	      if ( pk3 == 'y' ) printf("pb-dot, ");
	      if ( pk4 == 'y' ) printf("r, ");
	      if ( pk5 == 'y' ) printf("s, ");
	      if ( pk6 == 'y' ) printf("R, ");
	      printf(" and nothing else.\n ");
	      
	      /* If any variable is in prediction mode, act upon it */
	      
	      printf("\n Predicting distribution of the mass of the pulsar, of the companion, ");
	      if ( pk2 == 'p' )
		{
		  printf("gamma, ");
		  FGAMMA = fopen("pdf_gamma.dat","w");
		}
	      if ( pk3 == 'p' )
		{
		  printf("pb-dot, ");
		  FPBDOT = fopen("pdf_pbdot.dat","w");
		}
	      if ( pk4 == 'p' )
		{
		  printf("r, ");
		  FRS = fopen("pdf_r.dat","w");
		}
	      if ( pk5 == 'p' )
		{
		  printf("s, ");
		  FSS = fopen("pdf_ss.dat","w");
		}
	      if ( pk6 == 'p' )
		{
		  printf("R, ");
		  FSS = fopen("pdf_R.dat","w");
		}
	      printf(" and nothing else.\n ");
	      
	      if ( pk1 == 'y' )
		{
		  
		  /* start cycle of odots */
		  
		  odot = wdot - 5*ewdot;
		  dodot = ewdot/100.0;
		  
		  while (odot < ( wdot + 5*ewdot) )
		    {
		      
		      odot = odot + dodot;
		      
		      /* printf("%e \n", odot); */
		      
		      M = odot * (1 - e * e) * pow ( n , -five3rd) * pow( TSUN , -two3rd) / 3;
		      M = pow (M,1.5);
		      
		      /* what is the probability of seeing this odot? */
		      height0 = exp (- (odot - wdot)*(odot - wdot)/(2*ewdot*ewdot));
		      
		      /* printf("\n M = %f, odot = %e, height = %f \n", M, odot, height); */
		      
		      /* Now: go through the sin i loop :
			 probabilities are constant in the cos i domain,
			 so that should be transformed into the sin i domain, as to increase
			 the spead of calculation  */
		      
		      for (sini = dsin; sini < 1.0000; sini = sini + 2*dsin)
			{
			  cosi = sqrt(1 - sini*sini);
			  m2 = pow (f, one3rd) * pow(M, two3rd) / sini;
			  m1 = M - m2;
		      
			  gamma = gammak*m2*(M+m2)*pow(M,-four3rd);
			  pbdot = pbdotk* m1 * m2 / pow(M, one3rd);
			  rs = m2;
			  ss = sini;
			  R = m1/m2;
			  
			  /* printf("\n cos i = %f, sin i = %f, m1 = %f, m2 = %f, \n", sini, sini, m1, m2); */
			  
			  /* Now: check the bin for m1 and m2 */
			  j = floor(NM * 1.0 * ((m1-m1l)/(m1u-m1l)));
			  k = floor(NM * 1.0 * ((m2-m2l)/(m2u-m2l)));
			  
			  /* DO THE DRILL: TAKE INTO ACCOUNT WHATEVER OTHER CONSTRAINTS HAVE BEEN SPECIFIED TO BE IMPORTANT */
			  
			  
			  if ( pk0 == 'y' )
			    {
			      /* Update numbers taking into account probability for inclination */
			      sin1 = sini - dsin;
			      sin2 = sini + dsin;
			      cos1 = sqrt(1 - sin1*sin1);
			      cos2 = sqrt(1 - sin2*sin2);

			      height = height0 * 10000.0 * (cos1 - cos2);
			    }
			  else
			    height = height0;
			  
			  if ( pk2 == 'y' )
			    {
			      
			      /* Update numbers taking into account gamma */
			      exponent = - ( gamma - gammam ) * ( gamma - gammam )/( 2 * egamma * egamma );
			      height = height * exp(exponent); 
			    }
			  
			  if ( pk3 == 'y' )
			    {
			      
			      /* Update numbers taking into account pbdot */
			      exponent = - ( pbdot - pbdotm ) * ( pbdot - pbdotm )/( 2 * epbdot * epbdot );
			      height = height * exp(exponent);
			      
			    }
			  
			  if ( pk4 == 'y' )
			    {
			      
			      /* Update numbers taking into account r */
			      exponent = - ( rs - rsm ) * ( rs - rsm )/( 2 * ers * ers );
			      height = height * exp(exponent);
			      
			    }
			  
			  if ( pk5 == 'y' )
			    {
			      
			      /* Update numbers taking into account s */
			      exponent = - ( ss - ssm ) * ( ss - ssm )/( 2 * ess * ess );
			      height = height * exp(exponent);
			  
			    }
			  
			  if ( pk6 == 'y' )
			    {
			      
			      /* Update numbers taking into account s */
			      exponent = - ( R - Rm ) * ( R - Rm )/( 2 * eR * eR );
			      height = height * exp(exponent);			  
			    }

			  /* AND THEN PROJECT THE pdfs */
			  
			  if ( (j > -1) && (k > -1) )
			    {
			      mtotal = mtotal + height;
			      if ((j < NM) && (k < NM) )
				{
				  i = j*NM + k;
				  npc[i] = npc[i] + height;
				  np[j] = np[j] + height;
				  nc[k] = nc[k] + height;
				  
				}
			      
			      if ( pk2 == 'p' )
				{
				  l = floor(NM * (gamma - gammam + 5* egamma) /(10 * egamma));
				  if ((l > -1) && (l < NM)) fgamma[l] = fgamma[l] + height;
				}
			      
			      if ( pk3 == 'p' )
				{
				  l = floor(NM * (pbdot - pbdotm + 5* epbdot) /(10 * epbdot));
				  if ((l > -1) && (l < NM)) fpbdot[l] = fpbdot[l] + height;
				}
			      
			      if ( pk4 == 'p' )
				{
				  l = floor(NM * (rs - rsm + 5* ers) /(10 * ers));
				  if ((l > -1) && (l < NM)) frs[l] = frs[l] + height;
				}
			      
			      if ( pk5 == 'p' )
				{
				  l = floor(NM * (ss - ssm + 5* ess) /(10 * ess));
				  if ((l > -1) && (l < NM)) fss[l] = fss[l] + height;
				}
			      
			      if ( pk6 == 'p' )
				{
				  l = floor(NM * (R - Rm + 5* eR) /(10 * eR));
				  if ((l > -1) && (l < NM)) fR[l] = fR[l] + height;
				}
			      
			    }
			}
		    }
		}
	      
	      /* end of the omega-dot loop */
	      
	      /* if we have no omega-dot, then we assume that the other constraints can be seen on the map - do a grid search */
	      
	      else
		{
		  
		  /* grid elements */

		  dm1 = (m1u - m1l)/(NM * 1.0);
		  dm2 = (m2u - m2l)/(NM * 1.0);
		  
		  /* calculate masses and sin i */
		  for (j = 0; j < NM; j++)
		    {
		      m1 = m1l + j * dm1;
		      
		      for (k = 0; k < NM; k++)
			{
			  /* calculate sin i */
			  
			  m2 = m2l + k * dm2;
			  
			  M = m1 + m2;
			  
			  sini = pow ((f * M * M), one3rd) / m2;
			  
			  /* calculate by how much the angle would change if dm is added */
			  sinidiff = pow ((f * (M+0.000001) * (M+0.000001)), one3rd) / (m2+0.000001);
			  
			  height = 1000000.0;
			  
			  /* if we are in the region allowed by the mass function, then calculate */
			  
			  if (sini < 1)
			    {
			      
			      cosi = sqrt(1 - sini*sini);
			      cosidiff = sqrt(1 - sinidiff*sinidiff);
			      odot = pow ( n , five3rd) / (1 - e * e) * pow( (M*TSUN) , two3rd) * 3;
			      gamma = gammak*m2*(M+m2)*pow(M,-four3rd);
			      pbdot = pbdotk* m1 * m2 / pow(M, one3rd);
			      rs = m2;
			      ss = sini;
			      R = m1/m2;
			      
			      /* DO THE DRILL AGAIN: TAKE INTO ACCOUNT WHATEVER OTHER CONSTRAINTS HAVE BEEN SPECIFIED TO BE IMPORTANT */
			      
			      if ( pk0 == 'y' )
				{
				  /* Update numbers taking into account probability for inclination */
				  height = height * (cosidiff - cosi);
				}
			      
			      /* if ( pk1 == 'y' )
				{
				  Update numbers taking into account probability for inclination
				  exponent = - ( odot - wdot ) * ( odot - wdot ) /( 2 * ewdot * ewdot );
				  height = height * exp(exponent);
				} */

			      if ( pk2 == 'y' )
				{
				  
				  /* Update numbers taking into account gamma */
				  exponent = - ( gamma - gammam ) * ( gamma - gammam )/( 2 * egamma * egamma );
				  height = height * exp(exponent);
				  
				}
			      
			      if ( pk3 == 'y' )
				{
				  
				  /* Update numbers taking into account pbdot */			      
				  exponent = - ( pbdot - pbdotm ) * ( pbdot - pbdotm )/( 2 * epbdot * epbdot );
				  height = height * exp(exponent);
				  
				}
			      
			      if ( pk4 == 'y' )
				{
				  
				  /* Update numbers taking into account r */
				  exponent = - ( rs - rsm ) * ( rs - rsm )/( 2 * ers * ers );
				  height = height * exp(exponent);
				  
				}
			      
			      if ( pk5 == 'y' )
				{
				  
				  /* Update numbers taking into account s */
				  exponent = - ( ss - ssm ) * ( ss - ssm )/( 2 * ess * ess );
				  height = height * exp(exponent);
				  
				}
			      
			      if ( pk6 == 'y' )
				{
				  
				  /* Update numbers taking into account s */
				  exponent = - ( R - Rm ) * ( R - Rm )/( 2 * eR * eR );
				  height = height * exp(exponent);			  
				}
			      
			      
			      /* AND THEN PROJECT THE pdfs */
			      
			      if ( (j > -1) && (k > -1) )
				{
				  mtotal = mtotal + height;
				  if ((j < NM) && (k < NM) )
				    {
				      i = j*NM + k;
				      npc[i] = npc[i] + height;
				      np[j] = np[j] + height;
				      nc[k] = nc[k] + height;
				    }
				  
				  if ( pk2 == 'p' )
				    {
				      l = floor(NM * (gamma - gammam + 5* egamma) /(10 * egamma));
				      if ((l > -1) && (l < NM)) fgamma[l] = fgamma[l] + height;
				    }
				  
				  if ( pk3 == 'p' )
				    {
				      l = floor(NM * (pbdot - pbdotm + 5* epbdot) /(10 * epbdot));
				      if ((l > -1) && (l < NM)) fpbdot[l] = fpbdot[l] + height;
				    }
				  
				  if ( pk4 == 'p' )
				    {
				      l = floor(NM * (rs - rsm + 5* ers) /(10 * ers));
				      if ((l > -1) && (l < NM)) frs[l] = frs[l] + height;
				    }
				  
				  if ( pk5 == 'p' )
				    {
				      l = floor(NM * (ss - ssm + 5* ess) /(10 * ess));
				      if ((l > -1) && (l < NM)) fss[l] = fss[l] + height;
				    }
				  
				  if ( pk6 == 'p' )
				    {
				      l = floor(NM * (R - Rm + 5* eR) /(10 * eR));
				      if ((l > -1) && (l < NM)) fR[l] = fR[l] + height;
				    }

				  
				}
			    } 		      
			}
		    }
		}
	    }
	  
	  /* The 2-d p.d.f. has been calculated! */
	  
	  /* first: let's print the predictions */
	  
	  if ( pk2 == 'p' )
	    {
	      mptotal = 0;
	      for (l = 0; l < NM; l++)
		{
		  gamma = gammam + (l * 1.0 - NM * 0.5) / (0.1 * NM) * egamma;
		  fprintf(FGAMMA," %f  %f  %f \n", gamma, fgamma[l], mptotal/mtotal);
		  mptotal = mptotal + fgamma[l];
		}
	      fclose(FGAMMA);
	    }
	  
	  if ( pk3 == 'p' )
	    {
	      mptotal = 0;
	      for (l = 0; l < NM; l++)
		{
		  pbdot = pbdotm + (l * 1.0 - NM * 0.5) / (0.1 * NM) * epbdot;
		  fprintf(FPBDOT," %e  %f  %f \n", pbdot, fpbdot[l], mptotal/mtotal);
		  mptotal = mptotal + fpbdot[l];
		}
	      fclose(FPBDOT);
	    }
	  
	  if ( pk4 == 'p' )
	    {
	      mptotal = 0;
	      for (l = 0; l < NM; l++)
		{
		  rs = rsm + (l * 1.0 - NM * 0.5) / (0.1 * NM) * ers;
		  fprintf(FRS," %f  %f  %f \n", rs, frs[l], mptotal/mtotal);
		  mptotal = mptotal + frs[l];
		}
	      fclose(FRS);
	    }
	  
	  if ( pk5 == 'p' )
	    {
	      mptotal = 0;
	      for (l = 0; l < NM; l++)
		{
		  ss = ssm + (l * 1.0 - NM * 0.5) / (0.1 * NM) * ess;
		  fprintf(FSS," %f  %f  %f \n", ss, fss[l], mptotal/mtotal);
		  mptotal = mptotal + fss[l];
		}
	      fclose(FSS);
	    }
	  
	  if ( pk6 == 'p' )
	    {
	      mptotal = 0;
	      for (l = 0; l < NM; l++)
		{
		  R = Rm + (l * 1.0 - NM * 0.5) / (0.1 * NM) * eR;
		  fprintf(FR," %f  %f  %f \n", R, fR[l], mptotal/mtotal);
		  mptotal = mptotal + fR[l];
		}
	      fclose(FR);
	    };
	  
	  /* do we want contour plots? */
	  
	  if (d2 == 'y')
	    {
	      
	      /* define targets for probability inside each contour */
	      
	      /* First contour: should include 99.74% of probability */
	      target[0] = 99.74; alev[0] = 50000000;
	      
	      /* Second contour should include 95.44% of probability */
	      target[1] = 95.44; alev[1] = 50000000;
	      
	      /* Third contour should include 68.26% of probability */
	      target[2] = 68.26; alev[2] = 50000000;
	      
	      /* let us define the hights of each contour */
	      
	      for (j = 0; j < NCONT; j++)
		{
		  c = 0;
		  k = 0;
		  
		  /* start contour from zero */
		  alevc[j] = 0;
		  
		  while ( c == 0 )
		    {
		      /* add all points above countour level */
		      for (i =0; i < NM2; i++)
			if (npc[i] > alev[j])
			  alevc[j] = alevc[j] + npc[i];
		      
		      /* transform this into a percentage */
		      alevc[j] = 100 * alevc[j]/mtotal;
		      
		      /* is their sum the percentage we want? */
		      if ( (fabs(alevc[j] - target[j])) < 0.005 )
			{
			  c = 1;
			}
		      else
			{
			  /* if not enough probability inside contour, lower contour level,
			     otherwise raise contour level */
			  if ( (alevc[j] - target[j]) < 0 )
			    alev[j] = alev[j] * 0.99;
			  else
			    alev[j] = alev[j] * 1.0015;
			  k++;
			  /* if we enter a wicked loop, stop */
			  if (k > 2000) c = 1;
			}
		    }
		}
	      
	      /* printf("\n Contour levels defined.\n");*/
	      
	      /* Now, the grand task: the countour plot */
	      
	      /* define transformation matrix */
	      tr[0] = m1l;
	      tr[1] = 0;
	      tr[2] = c1;
	      tr[3] = m2l;
	      tr[4] = c2;
	      tr[5] = 0;
	  
	      /* make floats out of the npc vector, for display purposes */
	      
	      for ( i = 0 ; i < NM2 ; i++ ) fnpc[i] = (float) npc[i];
	      
	      /* make contour plot */

	      cpgsls(1);
	      cpgslw(1);

	      for ( i = 0 ; i < NCONT ; i++ )
		{
		  cpgcont(fnpc,NM,NM,1,NM,1,NM,&alev[i],-1,tr);
		  printf("\n Contour number %d encloses %f percent of simulated points.\n", i, alevc[i]);
		}
	    }
	  
	  /* Do we want the 1-D p.d.f.s? */
	  
	  /* We are putting this after the contour plot so that the contour
	     plot falls in the right window! */
	  
	  if (d3 == 'y')
	    {
	      
	      /* If so, we should normalize p.d.f.s. */
	      
	      /* first: find maximum for probability density */
	      
	      npmax = 0; 
	      ncmax = 0;
	      
	      for (i = 0; i < NM; i++)
		{
		  if (np[i] > npmax) npmax = np[i];
		  if (nc[i] > ncmax) ncmax = nc[i];
		}
	      
	      /* specify targets for mass distribution */
	      
	      mtarget[0] = 4.56;
	      mtarget[1] = 50.00;
	      mtarget[2] = 95.44;
	  
	      mptotal = 0;
	      mctotal = 0;
      
	      /* open files where p.d.f.s are going to be written */
      
	      PM = fopen("pdf_m1.dat","w");
	      CM = fopen("pdf_m2.dat","w");
	      
	      /* find at which masses these targets are filled */
	      /* use same cycle to normalize and print */
	      
	      for (i = 0; i< NM; i++)
		{
		  
		  for (j = 0; j < NCONT; j++)
		    {
		      /* is total less than target ? */
		      flagp[j] = ((100.0*mptotal/mtotal) - mtarget[j]);
		      flagc[j] = ((100.0*mctotal/mtotal) - mtarget[j]);
		    }
		  
		  /* calculate cummulative numbers */
		  mptotal = mptotal + np[i];
		  mctotal = mctotal + nc[i];
		  
		  for (j = 0; j < NCONT; j++)
		    {
		      /* is total still less than target ? */
		      flagp[j] = flagp[j] * ((100.0*mptotal/mtotal) - mtarget[j]);
		      flagc[j] = flagc[j] * ((100.0*mctotal/mtotal) - mtarget[j]);
		  
		      if ( flagp[j] < 0 )
			{
			  /* if there was a change, let's record the number */
			  mpalev[j] = m1l + i * c1;
			  mpalevy[j] = np[i];
			  printf("\n Mass of pulsar for %f probability is %f\n", mtarget[j], mpalev[j]);
			}
		      
		      if ( flagc[j] < 0 )
			{
			  /* same for companion mass */
			  mcalev[j] = m2l + i * c2;
			  mcalevy[j] = nc[i];
			  printf("\n Mass of companion for %f probability is %f\n", mtarget[j], mcalev[j]);
			}
		    }
		  /* use opportunity to do linearization. Build here vector for plotting */
		  
		  fnp[i] = np[i]/npmax;
		  fnc[i] = nc[i]/ncmax;
		  fmp[i] = m1l + i * c1;
		  fmc[i] = m2l + i * c2;
	      
		  fprintf(PM," %f  %f  %f \n", fmp[i], fnp[i], mptotal/mtotal);
		  fprintf(CM," %f  %f  %f \n", fmc[i], fnc[i], mctotal/mtotal);
		}

	      fclose(PM);
	      fclose(CM);

	      /* before printing: normalize values for vertical lines in mass plots */
	      
	      for (j =0; j < NCONT; j++)
		{
		  mpalevy[j] = mpalevy[j]/npmax;
		  mcalevy[j] = mcalevy[j]/ncmax;
		}
	      
	      /* Open window to display distribution of pulsar mass */
	      
	      cpgsch(1.0);
	      cpgvsiz(1.0, 5.5, 5.6, 7.5);
	      cpgswin(m1l, m1u, 0, 1.1);
	      cpgbox("BCTS", 0.0, 0, "BCTS", 0.0, 0);
	      cpglab("", "Probability density", "");
	  
	      /* print pulsar mass distribution */
	      cpgline(NM, fmp, fnp);
	  
	      /* print lines for median */ 
	      mx[0] = mpalev[1];
	      mx[1] = mpalev[1];
	      my[0] = 0;
	      my[1] = mpalevy[1];
	      cpgline(2, mx, my);
	  
	      /* Open window to display distribution of companion mass */
	  
	      cpgvsiz(5.6, 7.5, 1.0, 5.5);
	      cpgswin(1.1, 0, m2l, m2u); 
	      cpgbox("BCTS", 0.0, 0, "BCTS", 0.0, 0);
	      cpglab("Probability density", "", "");
	  
	      /* print companion mass distribution */
	      cpgline(NM, fnc, fmc);
	      
	      /* print lines for median */ 
	      my[0] = mcalev[1];
	      my[1] = mcalev[1];
	      mx[0] = 0;
	      mx[1] = mcalevy[1];
	      cpgline(2, mx, my);
	      
	    }
	  
	  printf("\n");
      
	  /* close the figure */
	  cpgend();
      
	}
    }

  return EXIT_SUCCESS;

}