void GOMX3_Doses(){
  
  const int num_altitudes    = 5;
  const int num_inclinations = 4;
  const int num_points       = 70;

  ifstream Doses[num_altitudes][num_inclinations];  
  
  double Altitude[num_altitudes];
  double Inclination[num_inclinations];
  
  double TotalDose[num_altitudes][num_inclinations][num_points];  
  double Thickness[num_points];
  
  double Dose_1mm[num_altitudes][num_inclinations];
  double Dose_2mm[num_altitudes][num_inclinations];
  double Dose_5mm[num_altitudes][num_inclinations];
  
  TH2D *Hist_Dose1 = new TH2D("Dose_1","Total Dose (1mm Al) [krad];Orbit Altitude [km];Inclination [deg]",num_altitudes,400.,800.,num_inclinations,0.,90.);
  TH2D *Hist_Dose2 = new TH2D("Dose_2","Total Dose (2mm Al) [krad];Orbit Altitude [km];Inclination [deg]",num_altitudes,400.,800.,num_inclinations,0.,90.);
  TH2D *Hist_Dose5 = new TH2D("Dose_5","Total Dose (5mm Al) [krad];Orbit Altitude [km];Inclination [deg]",num_altitudes,400.,800.,num_inclinations,0.,90.);
  
  char File_Dose[20];
  int N;
  
  Altitude[0] = 400.;
  Altitude[1] = 500.;
  Altitude[2] = 600.;
  Altitude[3] = 700.;
  Altitude[4] = 800.;
  
  Inclination[0] = 0.;
  Inclination[1] = 30.;
  Inclination[2] = 60.;
  Inclination[3] = 90.;
  
  float A,B,C,D,E,F,G,H,I;
  
  for(int i=0;i<num_altitudes;i++){
    for(int j=0;j<num_inclinations;j++){
    
      N = sprintf(File_Dose,"%dH/%ddeg/dose.dat",Altitude[i],Inclination[j]);
      
      Doses[i][j].open(File_Dose);
      
      cout<<File_Dose<<endl;
      cout<<"Altitude = "<<Altitude[i]<<" km & Inclination = "<<Inclination[j]<<" deg"<<endl;
      //cout<<"Thickness(mm Al) Dose(rad)"<<endl;
      
      for(int k=0;k<num_points;k++){  
            //Xtrpe=0; Ztrpe=0;
        Doses[i][j]>>A>>B>>C>>D>>E>>F>>G>>H>>I;
        Thickness[k] = A;
        TotalDose[i][j][k] = I;    
        //cout<<k+1<<" "<<Thickness[k]<<"      "<<TotalDose[i][j][k]<<endl;
      }
      
      Dose_1mm[i][j] = 1e-3*TotalDose[i][j][10];
      Dose_2mm[i][j] = 1e-3*TotalDose[i][j][20];
      Dose_5mm[i][j] = 1e-3*TotalDose[i][j][32];
      
      Hist_Dose1->Fill(Altitude[i],Inclination[j],Dose_1mm[i][j]);
      Hist_Dose2->Fill(Altitude[i],Inclination[j],Dose_2mm[i][j]);
      Hist_Dose5->Fill(Altitude[i],Inclination[j],Dose_5mm[i][j]);
    }
  }
  
  TCanvas *c1 = new TCanvas("c1");
  c1->SetFillColor(10);
  c1->SetGrid();
  gStyle->SetOptStat(0);
  
  Hist_Dose1->Draw("cont4z");
      
  //legend = sprintf(buffer,"Radiance [W.m-2.sr-1.#mum-1] (Local Time=%.1f hours, #lambda=%.1f nm);Days Of the Year;Latitude [deg]",Hour1,chosen_Wavelength);     
  //TH2D *h = (TH2D *)gROOT->FindObject("Rad_31");
  //delete h;
          
}