// ========================================================================================//
// = Compilation: g++ -g `root-config --cflags` `root-config --libs` -o Reader Reader.cpp =//
// ========================================================================================//

#include <stdio.h>
#include <math.h>
#include <string>
#include <iostream>
#include <fstream>
#include "TRint.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TAxis.h"



//Class declaration:

class Map
{
public:
  int nb;
  TGraph* B;
  TGraph* k1;

  Map();
  Map(const char* file);
  ~Map();
};



//Main function;:

///int main(int argc, char *argv[])
void Reader()
{
  const char* filename = "Test_data.txt";
  
///  TApplication * app=new TApplication("ROOT Example",&argc, argv, NULL, 0);
  gStyle -> SetCanvasColor(0);
  gStyle -> SetCanvasBorderMode(0);
  gStyle -> SetFrameBorderMode(0);
  gStyle -> SetTitleSize(0.04);

  Map* Vmap = new Map(filename);

  TCanvas * can = new TCanvas("Vertical Map",filename,0,0, 1000, 1000);
  can -> Divide(1,2);
  can -> cd(1);
  Vmap -> B -> Draw("APL");
  can -> cd(2);
  gPad->SetLogy(); //Uncomment to see the problem.
  Vmap -> k1 -> Draw("APL");
  can->Modified();
  can->Update();
//  app->Run();  
//  return 0;
}

// ============================================================================
// ============================================================================

//Class used:


Map::Map() {}

Map::Map(const char* file)
{

  double bxbefore = 0;
  double bybefore = 0;
  double bzbefore = 0;
  double zbefore = 0;

  FILE *input;
  char skip[750];
  nb = -1;

 input = fopen(file, "r");

 while(feof(input)==false)
    {
      fgets(skip, 750, input);
      nb++;
    }
 
  fclose(input);

  B      = new TGraph(nb);
  k1     = new TGraph(nb-1);


  
  //Read the datas, do some calculations and set points:

  input = fopen(file, "r");
  
  char hour;
  double z, bx, by, bz;
  
  
  for(int n=0; n<nb; n++)
    {
      fscanf(input,"%s %lf %lf %lf %lf", 
	     &hour, &z, &bx, &by, &bz);    
      
      double modB = sqrt(bx*bx + by*by + bz*bz);
      
      double f = (bx*(bxbefore - bx) + by*(bybefore - by) + bz*(bzbefore - bz))/modB;
      double den1 = (bxbefore - bx - bx*f/modB)*(bxbefore - bx - bx*f/modB);
      den1 += (bybefore - by - by*f/modB)*(bybefore - by - by*f/modB);
      den1 += (bzbefore - bz - bz*f/modB)*(bzbefore - bz - bz*f/modB);
      den1 = sqrt(den1);
      double den2 = acos((bx*bxbefore + by*bybefore + bz*bzbefore)/modB/sqrt(bxbefore*bxbefore+bybefore*bybefore+bzbefore*bzbefore));
      double k1_calc=2.621e-4*modB*modB*(z-zbefore)/den1;
      double k2_calc=2.621e-4*modB*(z-zbefore)/den2;

      B      -> SetPoint(n, z, modB);
 
      if(n>0)
	{
	  k1     -> SetPoint(n-1, z, k1_calc);
	  
	}
      
      bxbefore = bx;
      bybefore = by;
      bzbefore = bz;
      zbefore  = z;
      
      
    }
  
  fclose(input);

  //Graphs properties (where the problem is for k1 graph)

  B         -> SetTitle("");
  k1        -> SetTitle("");

  B         -> SetMarkerStyle(20);
  k1        -> SetMarkerStyle(20);

  B         -> GetXaxis() -> SetTitle("X axis");
  k1        -> GetXaxis() -> SetTitle("X axis");

  B         -> GetXaxis() -> SetTitleSize(0.05);
  k1        -> GetXaxis() -> SetTitleSize(0.05);
 
  B         -> GetYaxis() -> SetTitle("First value");
  k1        -> GetYaxis() -> SetTitle("Second value");

  B         -> GetYaxis() -> SetTitleSize(0.05);
  k1        -> GetYaxis() -> SetTitleSize(0.05);

  B         -> GetXaxis() -> SetLabelSize(0.05);
  k1        -> GetXaxis() -> SetLabelSize(0.05);

  B         -> GetYaxis() -> SetLabelSize(0.05);
  k1        -> GetYaxis() -> SetLabelSize(0.05);

}