Index: trunk/Mars/hawc/atmabs.C
===================================================================
--- trunk/Mars/hawc/atmabs.C	(revision 19764)
+++ trunk/Mars/hawc/atmabs.C	(revision 19764)
@@ -0,0 +1,174 @@
+void atmabs()
+{
+    // =============================================================================
+    // measured flux (lambda) = emitted flux (lambda) * trans(H, lambda)
+    // =============================================================================
+    TH2F hist1 ("TransCeres",         "Ceres Transmission",            700-250, 250, 700, 50, 0, 50);
+    TH2F hist1a("TransCeresRayleigh", "Ceres Transmission (Rayleigh)", 700-250, 250, 700, 50, 0, 50);
+    TH2F hist1b("TransCeresOzone",    "Ceres Transmission (Ozone)",    700-250, 250, 700, 50, 0, 50);
+    TH2F hist1c("TransCeresAerosol",  "Ceres Transmission (Aerosol)",  700-250, 250, 700, 50, 0, 50);
+
+    const char *ozone    = "resmc/atmosphere-ozone.txt";
+    const char *aerosols = "resmc/atmosphere-aerosols.txt";
+
+    MAtmosphere atm(1, ozone, aerosols);
+
+    cout << atm.HasValidAerosol() << " " << atm.HasValidOzone() << endl;
+
+    for (int h=0; h<50; h++)
+    {
+        for (int l=250; l<700; l++)
+        {
+            MPhotonData ph;
+            ph.SetWavelength(l);
+            ph.SetProductionHeight(h*1e3*1e2);
+
+            hist1a.Fill(l, h, atm.CalcTransmission(h*1e5, l, 0));
+            hist1b.Fill(l, h, atm.CalcOzoneAbsorption(h*1e5, l, 0));
+            hist1c.Fill(l, h, atm.CalcAerosolAbsorption(h*1e5, l, 0));
+            hist1.Fill(l, h, atm.GetTransmission(ph));
+        }
+    }
+
+    // =============================================================================
+    // The CCD Photometric Calibration Cookbook
+    // 8 Atmospheric Extinction and Air Mass
+    // http://star-www.rl.ac.uk/star/docs/sc6.htx/sc6se8.html
+    //
+    // measured magnitude (lambda) = emitted magnitude (lambda) + extco(lambda) * X(z)
+    //
+    //  ==> measured flux = emitted flux * exp(extco * X(h))
+    //
+    // =============================================================================
+
+    TH2F hist2("AtmAbs", "Corsika Atmosphere atmabs.txt", (700-180)/5+1, 178.5, 702.5, 51, -0.5, 50.5);
+
+    ifstream fin("/home/tbretz/SW/corsika-76900/run/atmabs.dat");
+
+    TString buf;
+    buf.ReadLine(fin);
+    cout << buf << endl;
+
+    while (1)
+    {
+        double wl;
+        fin >> wl;
+        if (!fin)
+            break;
+
+        for (int h=0; h<51; h++)
+        {
+            double cc;
+            fin >> cc;
+            hist2.Fill(wl, h, cc);
+        }
+    }
+
+    // =============================================================================
+    hist1.SetMinimum(0);
+    hist1.SetMaximum(1);
+    hist1.SetContour(500);
+
+    hist1a.SetMinimum(0);
+    hist1a.SetMaximum(1);
+    hist1a.SetContour(500);
+
+    hist1b.SetMinimum(0);
+    hist1b.SetMaximum(1);
+    hist1b.SetContour(500);
+
+    hist1c.SetMinimum(0);
+    hist1c.SetMaximum(1);
+    hist1c.SetContour(500);
+
+    // =============================================================================
+
+    gStyle->SetOptStat(11);
+
+    TCanvas *c = new TCanvas;
+    c->Divide(2, 2);
+
+    c->cd(1);
+    hist1a.DrawCopy("colz");
+
+    c->cd(2);
+    hist1b.DrawCopy("colz");
+
+    c->cd(3);
+    hist1c.DrawCopy("colz");
+
+    c->cd(4);
+    hist1.DrawCopy("colz");
+
+    // =============================================================================
+
+    c = new TCanvas;
+    c->Divide(2, 2);
+
+    c->cd(1);
+    hist1.DrawCopy("colz");
+
+    c->cd(2);
+    hist2.SetMinimum(0);
+    //hist2.SetMaximum(1);
+    hist2.SetContour(500);
+    hist2.DrawCopy("colz");
+
+    // =============================================================================
+    //
+    //  see CORSIKA ATABSO()
+    //
+    //  CORSIKA calculates a coefficient COATEX which is the value interpolated
+    //  in wavelength and height from the table.
+    //
+    //  From this value, the interpolated value at observation level is subtraced.
+    //
+    //  The probability to survive (=transmission) is then exp(-COATEX/WEMIS)
+    //
+    //  WEMIS defined as the particles z-direction cosine (-> correction for
+    // zenith angle (flat atmosphere))
+    //
+    // =============================================================================
+
+    TH2F hist4("TransCorsika", "Corsika Transmission", (700-180)/5+1, 178.5, 702.5, 51, -0.5, 50.5);
+    for (int x=1; x<=hist2.GetNbinsX(); x++)
+    {
+        double sum = 0;
+        for (int y=1; y<=hist2.GetNbinsY(); y++)
+        {
+            const double cosz = 1;
+            hist4.SetBinContent(x, y, exp(-hist2.GetBinContent(x, y)/cosz));
+        }
+    }
+
+    c->cd(4);
+    hist4.SetMinimum(0);
+    hist4.SetMaximum(1);
+    hist4.SetContour(500);
+    hist4.DrawCopy("colz");
+
+    c->cd(3);
+    TH1D *p1 = hist1.ProjectionX("TransCeres50km",   hist1.GetNbinsY(), hist1.GetNbinsY());
+    TH1D *p4 = hist4.ProjectionX("TransCorsika50km", hist4.GetNbinsY(), hist4.GetNbinsY());
+    TH1D *p1b = hist1.ProjectionX("TransCeres5km",   5, 5);
+    TH1D *p4b = hist4.ProjectionX("TransCorsika5km", 5, 5);
+
+    p1->SetNameTitle("Trans", "Transmission");
+
+    p4->SetLineColor(kBlue);
+    p4b->SetLineColor(kBlue);
+    p1->SetMinimum(0);
+    p1->SetMaximum(1);
+    p1->Draw();
+    p4->Draw("same");
+    p1b->Draw("same");
+    p4b->Draw("same");
+
+    TLatex text;
+    text.DrawLatex(275, 0.850, "Ceres");
+    text.DrawLatex(340, 0.250, "50km");
+    text.DrawLatex(320, 0.650, "5km");
+    text.SetTextColor(kBlue);
+    text.DrawLatex(275, 0.925, "Corsika");
+
+}