// Functions associated with map projections...
// -----------------------------------------------------------------------------
/// \brief Transcendental equation for intermediate Mollweide parameter gamma.
// -----------------------------------------------------------------------------
Double_t
f(Double_t g, Double_t th)
{
   return (2.0 * g + sin(2.0 * g) - TMath::Pi() * sin(th));
}

// -----------------------------------------------------------------------------
/// \brief First derivative of eq. for intermediate Mollweide parameter gamma.
// -----------------------------------------------------------------------------
Double_t
df(Double_t g, Double_t th)
{
   return (2.0 + 2.0 * cos(2.0 * g));
}

// -----------------------------------------------------------------------------
/// \brief Newton-Raphson root finding algorithm used to solve for gamma.
// -----------------------------------------------------------------------------
Double_t
root_find(
   Double_t (*f)(Double_t, Double_t),
   Double_t (*df)(Double_t, Double_t),
   Double_t a,
   Double_t b,
   Double_t th
   )
{
   static Double_t precision = 1.0e-4;
   Double_t dx = b - a;
   Double_t x0 = (a + b) / 2.0;
   Double_t x1;

   do {
      x1 = x0 - f(x0, th) / df(x0, th);
      dx = x1 - x0;
      x0 = x1;
   } while (fabs(dx) > precision);

   return x0;
}

// -----------------------------------------------------------------------------
/// \brief Mollweide projection: long,lat --> x,y
// -----------------------------------------------------------------------------
void
project(const Double_t l, const Double_t b, Double_t *x, Double_t *y)
{
  Double_t g = root_find(f, df, -0.5*TMath::Pi(), 0.5*TMath::Pi(), b);
  *x = 2. * TMath::Sqrt(2.) * l * TMath::Cos(g) / TMath::Pi();
  *y = TMath::Sqrt(2.) * TMath::Sin(g);
}

// -----------------------------------------------------------------------------
/// \brief "Main" program.
// -----------------------------------------------------------------------------
void
PlotCoords(const Char_t *file, const Int_t nlon=9, const Int_t nlat=7)
{
  ifstream ifstr(file);
  if (!ifstr) {
    cerr << "\n Could not open " << file << "\n\n";
    return;
  }

  const Double_t deg = TMath::DegToRad();

  Double_t dlat = 180.*deg / (nlat - 1);
  Double_t dlon = 360.*deg / (nlon - 1);

  Int_t count, line;
  Double_t l, b, x, y;

  TGraph *grlon   = new TGraph[nlon];
  TGraph *grlat   = new TGraph[nlat];
  TGraph *grdata  = new TGraph;
  TMultiGraph *mg = new TMultiGraph;

  // Import data from a file.
  count = 0;
  while (ifstr >> b >> l) {
                                          // Galactic coords: make sure l > 180
    l = l>180. ? -(l-360.) : -l;          // runs from -180 to 0, and the plot
                                          // goes +180 +90 0 -90 -180.
    project(l*deg, b*deg, &x, &y);
    grdata->SetPoint(count++, x, y);
  }
  mg->Add(grdata, "p");
  grdata->SetMarkerColor(9);
  grdata->SetMarkerStyle(4);
  grdata->SetMarkerSize(0.8);

  // Lines of longitude.
  line = 0;
  for (l = -180*deg; l <= 180*deg; l += dlon) {
    count = 0;
    for (Double_t sinb = -1.; sinb <= 1.; sinb += 0.2) {
      b = TMath::ASin(sinb);
      project(l, b, &x, &y);
      grlon[line].SetPoint(count++, x, y);
    }
    grlon[line].SetLineColor(13);
    mg->Add(&grlon[line], "c");
    ++line;
  }

  // Lines of latitude.
  line = 0;
  for (b = -90*deg; b <= 90*deg; b += dlat) {
    count = 0;
    for (l = -180*deg; l <= 180*deg; l += 12*deg) {
      project(l, b, &x, &y);
      grlat[line].SetPoint(count++, x, y);
    }
    grlat[line].SetLineColor(13);
    mg->Add(&grlat[line], "c");
    ++line;
  }

  // Draw everything with a nice aspect ratio.
  TCanvas *c = new TCanvas("c", "c", 10., 10., 950, 500);
  c->SetFrameLineColor(0);
  c->SetFrameBorderSize(0);
  c->SetLeftMargin(0.05);
  c->SetRightMargin(0.05);
  c->SetTopMargin(0.05);
  c->SetBottomMargin(0.05);
  
  gStyle->SetOptStat(0);
  TH1D *hframe = new TH1D("", "", 10, -2.83, 2.83);
  hframe->SetMinimum(-1.42);
  hframe->SetMaximum(+1.42);
  TAxis *axis = hframe->GetXaxis();
  axis->SetAxisColor(0);
  axis->SetLabelColor(0);
  axis = hframe->GetYaxis();
  axis->SetAxisColor(0);
  axis->SetLabelColor(0);
  hframe->Draw();

  TPaveText *pt = new TPaveText(-3.15, -0.1, -2.85, 0.1);
  pt->SetLineColor(0);
  pt->SetFillColor(0);
  pt->SetTextSize(0.04);
  pt->SetTextFont(42);
  pt->AddText("+180#circ");
  pt->SetBorderSize(0);
  pt->Draw();

  pt = new TPaveText(2.85, -0.1, 3.15, 0.1);
  pt->SetLineColor(0);
  pt->SetFillColor(0);
  pt->SetTextSize(0.04);
  pt->SetTextFont(42);
  pt->AddText("-180#circ");
  pt->SetBorderSize(0);
  pt->Draw();

  mg->Draw("c");
}