/* ======================================================================== *\
! $Name: not supported by cvs2svn $:$Id: plotoptical.C,v 1.6 2006-11-17 11:45:29 dorner Exp $
! --------------------------------------------------------------------------
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz, 05/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Daniela Dorner, 05/2005 <mailto:dorner@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2006
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// plotoptical.C
// =============
//
// This macro is used to read optical data from the DB and plot them.
// 
// In the DB these values are stored in the tables Calibration and Star.
//
//
// To plot the whole database simple use:
//   .x plotdb.C+
//
// Plot only data for one source. For the available sources see the database
//   .x plotdb.C+("source")
//
// You can chose are certain MAGIC-period:
//   .x plotdb.C+(25)   --> all values from period 25 are plotted
//
// MAGIC periods correspond to one moon-phase and are defined as
// moon period-284. For details see MAstro::MoonPeriod and
// MAstro::MagicPeriod.
//
// or a time period from a certain date to a certain date
//   .x plotdb.C+("2004-11-14 00:00:00", "2005-02-28 00:00:00")
//  --> all values from 14.11.2004 0h to 28.2.2005 0h are plotted
//   .x plotdb.C+("2004-11-14 00:00:00", "2005-02-28 00:00:00", "source")
//  --> all values from 14.11.2004 0h to 28.2.2005 0h of "source" are plotted
//
//
// For details of the plots produced see the function plotall() below.
//
//
// The plot title and axis-titles are created by:
//    plot.SetDescription("Title;x", "TabName");
//
// Drawing the plot is initiated by
//    plot.Plot("OpticalData.fExposure", min, max, width);
//
// While OpticalData.fExposure can be any kind of variable to plot.
// min and max are the minimum and maximum of the histogram which is
// filled and width is the bin-width of this histogram.
//
// To group data (average) of a certain period use:
//    plot.GroupBy(MPlot::kGroupByNight);
// before initiating the plot.
//
//
// Make sure, that database and password are corretly set in a resource
// file called sql.rc and the resource file is found.
//
/////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <iomanip>

#include <TH1.h>
#include <TEnv.h>
#include <TPad.h>
#include <TLine.h>
#include <TText.h>
#include <TFrame.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TPRegexp.h>
#include <TSQLRow.h>
#include <TSQLResult.h>
#include <TGraphErrors.h>

#include "MTime.h"
#include "MAstro.h"
#include "MDataSet.h"
#include "MSQLMagic.h"
#include "MStatusDisplay.h"

class MPlot : public MParContainer
{
public:
    // Possible constants to group-by (average) over a certain period
    enum GroupBy_t
    {
        kNone,
        kGroupByPrimary,
        kGroupByHour,
        kGroupByNight,
        kGroupByWeek,
        kGroupByMonth,
        kGroupBySeason,
        kGroupByYear
    };

private:
    MSQLMagic &fServer;         // Reference to the sql-server class

    MDataSet *fDataSet;         // A possible dtaset to highlite single points

    TString   fPrimaryDate;     // The name of the data we plot
    TString   fPrimaryNumber;   // The corresponding name for the key number
    TString   fSecondary;       // The value versus which the second plot is made

    TString   fRequestFrom;     // Start of a requested date range
    TString   fRequestTo;       // End of a requested date range
    Int_t     fRequestPeriod;   // A possible requested period

    Float_t   fPlotMin;
    Float_t   fPlotMax;

    Float_t   fHistMin;
    Float_t   fHistMax;

    TString   fDescription;     // The description (title) of the plot
    TString   fNameTab;         // The name of the tab in the display

    TString   fCondition;       // An additional condition added to the query
    GroupBy_t fGroupBy;         // A possible Group-By flag

    // --------------------------------------------------------------------------
    //
    // Function to plot the result of the query
    //
    void PlotTable(TSQLResult &res, TString name, Float_t fmin, Float_t fmax, Float_t resolution)
    {
        // Enable all otions in the statistics box
        gStyle->SetOptStat(111111);

        TSQLRow *row;

        // Create TGraph objects
        TGraph &gt = res.GetFieldCount()>4 ? *new TGraphErrors : *new TGraph;
        gt.SetNameTitle(name, Form("%s vs Time", name.Data()));
        gt.SetMarkerStyle(kFullDotMedium);

        TGraph gz;
        gz.SetNameTitle(name, Form("%s vs <Zd>", name.Data()));
        gz.SetMarkerStyle(kFullDotMedium);

        TGraph gt0, gt1;
        gt0.SetMarkerColor(kRed);
        gt1.SetMarkerColor(kBlue);
        gt0.SetMarkerStyle(kFullDotLarge);
        gt1.SetMarkerStyle(kFullDotLarge);

        TGraph gz0, gz1;
        gz0.SetMarkerColor(kRed);
        gz1.SetMarkerColor(kBlue);
        gz0.SetMarkerStyle(kFullDotLarge);
        gz1.SetMarkerStyle(kFullDotLarge);

        Int_t first = -1;
        Int_t last  = -1;

        // Loop over the data
        while ((row=res.Next()))
        {
            // Get all fields of this row
            const char *date = (*row)[0];
            const char *zd   = (*row)[1];
            const char *val  = (*row)[2];
            const char *snum = res.GetFieldCount()>3 ? (*row)[3] : 0;
            const char *verr = res.GetFieldCount()>4 ? (*row)[5] : 0;
            if (!date || !val || !zd)
                continue;

            // check if date is valid
            MTime t(date);
            if (!t.SetSqlDateTime(date))
                continue;

            // check if it belongs to the requested MAGIC period
            if (fRequestPeriod>0 && MAstro::GetMagicPeriod(t.GetMjd())!=fRequestPeriod)
                continue;

            // Get axis range
            if (first<0)
                first = TMath::Nint(TMath::Floor(t.GetMjd()));
            last = TMath::Nint(TMath::Ceil(t.GetMjd()));

            // Convert a possible key number into a integer
            UInt_t seq = snum ? atoi(snum) : 0;

            // convert primary and secondary value into floats
            Float_t value = atof(val);
            Float_t zenith = atof(zd);

            // If a datset is given add the point to the special TGraphs
            // used for highliting these dates
            if (fDataSet)
            {
                if (fDataSet->HasOnSequence(seq))
                {
                    gt1.SetPoint(gt1.GetN(), t.GetAxisTime(), value);
                    gz1.SetPoint(gz1.GetN(), zenith, value);
                }

                if (fDataSet->HasOffSequence(seq))
                {
                    gt0.SetPoint(gt0.GetN(), t.GetAxisTime(), value);
                    gz0.SetPoint(gz0.GetN(), zenith, value);
                }
            }

            // Add Data to TGraph
            gt.SetPoint(gt.GetN(), t.GetAxisTime(), value);
            gz.SetPoint(gz.GetN(), zenith, value);

            // Set error-bar, if one
            if (verr)
                static_cast<TGraphErrors&>(gt).SetPointError(gt.GetN()-1, 0, atof(verr));
        }

        // If this is done earlier the plots remain empty since root 5.12/00
        if (fmax>fmin)
        {
            gt.SetMinimum(fmin);
            gt.SetMaximum(fmax);
            gz.SetMinimum(fmin);
            gz.SetMaximum(fmax);
        }

        gROOT->SetSelectedPad(0);

        // Create a TCanvas or open a new tab
        TString title = fNameTab.IsNull() ? name(name.First('.')+2, name.Length()) : fNameTab;
        TCanvas &c = fDisplay ? fDisplay->AddTab(title) : *new TCanvas;
        // Set fillcolor, remove border and divide pad
        c.SetFillColor(kWhite);
        c.SetBorderMode(0);
        c.Divide(1,2);

        // Output mean and rms to console
        cerr << setprecision(4) << setw(10) << title << ":   ";
        if (gt.GetN()==0)
        {
            cerr << "     <empty>" << endl;
            return;
        }
        cerr << setw(8) << gt.GetMean(2) << "+-" << setw(8) << gt.GetRMS(2) << "   ";
        if (gt0.GetN()>0 || gt1.GetN()>0)
        {
            cerr << setw(8) << gt1.GetMean(2) << "+-" << setw(8) << gt1.GetRMS(2) << "   ";
            cerr << setw(8) << gt0.GetMean(2) << "+-" << setw(8) << gt0.GetRMS(2);
        }
        cerr << endl;

        TVirtualPad *pad = gPad;

        // draw contants of pad 2 (counting starts at 0)
        pad->cd(2);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        gPad->SetGridy();

        gPad->SetLeftMargin(0.06);
        gPad->SetRightMargin(0.06);
        gPad->SetBottomMargin(0.08);

        // format axis
        TH1 *h = gt.GetHistogram();

        h->SetXTitle("Time");
        h->SetYTitle(name);
        h->GetXaxis()->SetTimeDisplay(1);
        h->GetYaxis()->SetTitleOffset(0.8);
        h->GetXaxis()->SetTitleOffset(1.0);
        h->GetXaxis()->SetLabelOffset(0.01);

        // draw TGraph
        gt.DrawClone("AP");
        if (gt0.GetN()>0)
            gt0.DrawClone("P");
        if (gt1.GetN()>0)
            gt1.DrawClone("P");

        // Add lines and text showing the MAGIC periods
        TLine l;
        TText t;
        Int_t num=0;
        l.SetLineStyle(kDotted);
        l.SetLineColor(kBlue);
        t.SetTextColor(kBlue);
        l.SetLineWidth(1);
        t.SetTextSize(h->GetXaxis()->GetLabelSize());
        t.SetTextAlign(21);
        Int_t p0 = MAstro::GetMagicPeriod(first);
        for (Int_t p = first; p<last; p++)
        {
            Int_t p1 = MAstro::GetMagicPeriod(p);
            if (p1!=p0)
            {
                l.DrawLine(MTime(p).GetAxisTime(), h->GetMinimum(), MTime(p).GetAxisTime(), h->GetMaximum());
                t.DrawText(MTime(p+15).GetAxisTime(), h->GetMaximum(), Form("%d", p1));
                num++;
            }
            p0 = p1;
        }
        if (num<4)
            gPad->SetGridx();

        const Double_t min = fHistMin>fHistMax ? h->GetMinimum()-resolution/2 : fHistMin;
        const Double_t max = fHistMin>fHistMax ? h->GetMaximum()+resolution/2 : fHistMax;

        // Use this to save the pad with the time development to a file
        //gPad->SaveAs(Form("plotdb-%s.eps", title.Data()));

        // Go back to first (upper) pad, format it and divide it again
        pad->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        gPad->Divide(2,1);

        TVirtualPad *pad2 = gPad;

        // format left pad
        pad2->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();

        // Create histogram
        const Int_t n = resolution>0 ? TMath::Nint((max-min)/resolution) : 50;

        TH1F hist("Hist", Form("Distribution of %s", fDescription.IsNull() ? name.Data() : fDescription.Data()), n, min, max);
        hist.SetDirectory(0);

        // Fill data into histogra,
        for (int i=0; i<gt.GetN(); i++)
            hist.Fill(gt.GetY()[i]);

        // Format histogram
        if (fDescription.IsNull())
            hist.SetXTitle(name);
        hist.SetYTitle("Counts");

        // plot histogram
        hist.DrawCopy("");

        // format right pad
        pad2->cd(2);
        gPad->SetBorderMode(0);
        gPad->SetFrameBorderMode(0);
        gPad->SetGridy();

        // format graph
        TH1 *h2 = gz.GetHistogram();

        h2->SetXTitle("Zd");
        h2->SetYTitle(name);

        // draw graph
        gz.DrawClone("AP");

        if (gz0.GetN()>0)
            gz0.DrawClone("P");
        if (gz1.GetN()>0)
            gz1.DrawClone("P");
    }

public:
    MPlot(MSQLMagic &server) : fServer(server), fDataSet(NULL),
        fRequestPeriod(-1), fPlotMin(0), fPlotMax(-1), fHistMin(0), fHistMax(-1), fGroupBy(kNone)
    {
    }
    ~MPlot()
    {
        if (fDataSet)
            delete fDataSet;
    }
    void SetDataSet(const TString filename)
    {
        if (fDataSet)
        {
            delete fDataSet;
            fDataSet = NULL;
        }
        if (!filename.IsNull())
            fDataSet = new MDataSet(filename);
    }
    void SetPlotRange(Float_t min, Float_t max, Int_t n=5) { fPlotMin = min; fPlotMax = max; }
    void SetHistRange(Float_t min, Float_t max) { fHistMin = min; fHistMax = max; }
    void SetRequestRange(const char *from="", const char *to="") { fRequestFrom = from; fRequestTo = to; }
    void SetRequestPeriod(Int_t n=-1) { fRequestPeriod = n; }
    void SetCondition(const char *cond="") { fCondition = cond; }
    void SetDescription(const char *d, const char *t=0) { fDescription = d; fNameTab = t; }
    void SetGroupBy(GroupBy_t b=kGroupByWeek) { fGroupBy=b; }
    void SetPrimaryDate(const char *ts) { fPrimaryDate=ts; }
    void SetPrimaryNumber(const char *ts) { fPrimaryNumber=ts; }
    void SetSecondary(const char *ts) { fSecondary=ts; }

    Int_t QueryKeyOfSource(TString src)
    {
        return fServer.QueryKeyOfName("Object", src, kFALSE);
    }

    Bool_t Plot(const char *value, Float_t min=0, Float_t max=-1, Float_t resolution=0)
    {
        TString named  = fPrimaryDate;
        TString named2 = fSecondary;
        TString namev  = value;

        TString tablev = namev(0, namev.First('.'));
        TString valuev = namev(namev.First('.')+1, namev.Length());

        TString tabled = named(0, named.First('.'));
        TString valued = named(named.First('.')+1, named.Length());

        TString query="SELECT ";
        switch (fGroupBy)
        {
        case kNone:
        case kGroupByPrimary:
            query += valued;
            break;
        case kGroupByHour:
            query += Form("DATE_FORMAT(%s, '%%Y-%%m-%%d %%H:30:00') AS %s ", fPrimaryDate.Data(), valued.Data());
            break;
        case kGroupByNight:
            query += Form("DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y-%%m-%%d 00:00:00') AS %s ", fPrimaryDate.Data(), valued.Data());
            break;
        case kGroupByWeek:
            query += Form("DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%x%%v') AS %s ", fPrimaryDate.Data(), valued.Data());
            break;
        case kGroupByMonth:
            query += Form("DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y-%%m-15 00:00:00') AS %s ", fPrimaryDate.Data(), valued.Data());
            break;
        case kGroupBySeason:
            //query += Form("DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y-%%m-15 00:00:00') AS %s ", fPrimaryDate.Data(), valued.Data());
            break;
        case kGroupByYear:
            query += Form("DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y-08-15 00:00:00') AS %s ", fPrimaryDate.Data(), valued.Data());
            break;
        }

        if (fGroupBy==kNone)
        {
            query += ", ";
            query += fSecondary;
            query += ", ";
            query += value;
            query += ", ";
            query += fPrimaryNumber;
            query += " ";
        }
        else
        {
            query += ", AVG(";
            query += fSecondary;
            query += "), AVG(";
            query += value;
            query += "), ";
            query += fPrimaryNumber;
            query += ", STD(";
            query += fSecondary;
            query += "), STD(";
            query += value;
            query += ") ";
        }

        query += Form("FROM %s ", tabled.Data());

        TString where(fCondition);

        const Bool_t interval = !fRequestFrom.IsNull() && !fRequestTo.IsNull();
        if (interval)
        {
            if (!where.IsNull())
                where += " AND ";
            where += Form("%s BETWEEN '%s' AND '%s' ",
                          fPrimaryDate.Data(), fRequestFrom.Data(), fRequestTo.Data());
        }

        if (!where.IsNull())
            where += " AND ";
        where += fCondition;
        where += " ";

        // ------------------------------

        query += fServer.GetJoins(tabled, query+" "+where);

        if (!where.IsNull())
        {
            query += "WHERE ";
            query += where;
        }

        if (fGroupBy!=kNone)
        {
            query += Form("GROUP BY %s ", valued.Data());
            //query += Form(" HAVING COUNT(%s)=(COUNT(*)+1)/2 ", valuev.Data());
        }
        query += Form("ORDER BY %s ", valued.Data());

        // ------------------------------

        TSQLResult *res = fServer.Query(query);
        if (!res)
            return kFALSE;

        if (max>min)
            PlotTable(*res, namev, min, max, resolution);
        else
            PlotTable(*res, namev, fPlotMin, fPlotMax, resolution);

        delete res;
        return kTRUE;
    }
};

void plotall(MPlot &plot, TString source)
{
    // Setup here the values for timestamp and secondary (upper/right) plot
    plot.SetPrimaryDate("OpticalData.fTimestamp");
    plot.SetPrimaryNumber("OpticalData.fTimestamp");
    plot.SetSecondary("OpticalData.fZenithDistance");

    // This is the condition to take only the "ok" flagged data
    // and to restrict the query to a given source (if any)
    TString cond = "fStatusKEY=1";
    if (!source.IsNull())
    {
        const Int_t key = plot.QueryKeyOfSource(source);
        if (key<0)
            return;
        cond += Form(" AND Object.fObjectKEY=%d", key);

    }
    plot.SetCondition(cond);

    // Plot exposure
    plot.SetDescription("Exposure;T_{E} [s]", "Expo");
    plot.Plot("OpticalData.fExposure",  0, 900, 60);

    // plot sky level
    plot.SetDescription("Sky Level;B [s^{-1}]", "SkyLvl");
    plot.Plot("OpticalData.fSkyLevel/OpticalData.fExposure", 0, 5.0, 0.01);

    // plot FWHM
    plot.SetDescription("Full Width Half Maximum;FWHM [s^{-1}]", "Fwhm");
    plot.Plot("OpticalData.fFWHM/OpticalData.fExposure",  0, 0.05, 0.001);

    // plot Aperture Radius
    plot.SetDescription("Aperture Radius;R_{A}", "ApRad");
    plot.Plot("OpticalData.fApertureRadius", 0, 10, 1);

    /*
     // Plot instrumental magnitude
     plot.SetDescription("Instrumental Magnitude;M_{I}", "InstMag");
     plot.Plot("OpticalData.fInstrumentalMag",  0, 30, 0.5);

     // Plot error of instrumental magnitude
     plot.SetDescription("Instrumental Magnitude Error;\\sigma_{M}", "MagErr");
     plot.Plot("OpticalData.fInstrumentalMagErr",  0, 1, 0.01);
     */

    // Plot magnitude corrected for the exposure
    plot.SetDescription("m_{1};m_{1}", "M1");
    plot.Plot("OpticalData.fInstrumentalMag+2.5*log10(OpticalData.fExposure)", 10, 35, 0.2);

    // Now take out all points named */BL from further queries
    // And skip all sources the magnitude is not known
    cond += " AND Object.fObjectName NOT LIKE '%/BL' AND NOT ISNULL(Object.fMagnitude) ";
    plot.SetCondition(cond);

    // Formula to calculate the extinction
    TString ext("3080/25.0*pow(10, (OpticalData.fInstrumentalMag+2.5*log10(OpticalData.fExposure)-Object.fMagnitude)/-2.5)");
    // Add this to correct for the ZA dependancy
    //    ext += "+0.0028*fZenithDistance-0.08";

    // Group all data of one image together and plot extinction
    plot.SetGroupBy(MPlot::kGroupByPrimary);
    plot.SetDescription("m_{1}-m_{true} (Extinction per Image);m_{1}-m_{true}", "ExtImg");
    plot.Plot(ext, 0.05, 1.2, 0.01);

    // Group data hourly together and plot extinction
    plot.SetGroupBy(MPlot::kGroupByHour);
    plot.SetDescription("m_{1}-m_{true} (Extinction per Hour);m_{1}-m_{true}", "ExtHour");
    plot.Plot(ext, 0.5, 1.2, 0.01);

    // Group data hourly together and plot extinction
    plot.SetGroupBy(MPlot::kGroupByNight);
    plot.SetDescription("m_{1}-m_{true} (Extinction per Night);m_{1}-m_{true}", "ExtNight");
    plot.Plot(ext, 0.5, 1.2, 0.01);

    // Group data monthly together and plot extinction
    plot.SetGroupBy(MPlot::kGroupByMonth);
    plot.SetDescription("m_{1}-m_{true} (Extinction per Month);m_{1}-m_{true}", "ExtMonth");
    plot.Plot(ext, 0.5, 1.2, 0.01);

    // Group data yearly together and plot extinction
    plot.SetGroupBy(MPlot::kGroupByYear);
    plot.SetDescription("m_{1}-m_{true} (Extinction per Year);m_{1}-m_{true}", "ExtYear");
    plot.Plot(ext, 0.5, 1.2, 0.01);
}

int plotoptical(TString from, TString to, const char *source=0)
{
    TEnv env("sql.rc");

    MSQLMagic serv(env);
    if (!serv.IsConnected())
    {
        cout << "ERROR - Connection to database failed." << endl;
        return 0;
    }

    cout << "plotoptical" << endl;
    cout << "-----------" << endl;
    cout << endl;
    cout << "Connected to " << serv.GetName() << endl;
    cout << endl;

    MStatusDisplay *d = new MStatusDisplay;
    d->SetWindowName(serv.GetName());
    d->SetTitle(serv.GetName());

    MPlot plot(serv);
    // plot.SetDataSet(dataset);
    plot.SetDisplay(d);
    plot.SetRequestRange(from, to);
    plotall(plot, source);
    // Use this to create output plots automatically
    //  d->SaveAsRoot("plotoptical.root");
    //  d->SaveAsPS("plotoptical.ps");

    return 1;
}

int plotoptical(const char *source)
{
    TEnv env("sql.rc");

    MSQLMagic serv(env);
    if (!serv.IsConnected())
    {
        cout << "ERROR - Connection to database failed." << endl;
        return 0;
    }

    cout << "plotoptical" << endl;
    cout << "-----------" << endl;
    cout << endl;
    cout << "Connected to " << serv.GetName() << endl;
    cout << endl;

    MStatusDisplay *d = new MStatusDisplay;
    d->SetWindowName(serv.GetName());
    d->SetTitle(serv.GetName());

    MPlot plot(serv);
    // plot.SetDataSet(ds);
    plot.SetDisplay(d);
    plot.SetRequestRange("", "");
    plotall(plot, source);
    // Use this to create output plots automatically
    //  d->SaveAsRoot("plotoptical.root");
    //  d->SaveAsPS("plotoptical.ps");

    return 1;
}

int plotoptical(Int_t period, const char *source="")
{
    TEnv env("sql.rc");

    MSQLMagic serv(env);
    if (!serv.IsConnected())
    {
        cout << "ERROR - Connection to database failed." << endl;
        return 0;
    }

    cout << "plotoptical" << endl;
    cout << "-----------" << endl;
    cout << endl;
    cout << "Connected to " << serv.GetName() << endl;
    cout << endl;

    MStatusDisplay *d = new MStatusDisplay;
    d->SetWindowName(serv.GetName());
    d->SetTitle(serv.GetName());

    MPlot plot(serv);
    // plot.SetDataSet(dataset);
    plot.SetDisplay(d);
    plot.SetRequestPeriod(period);
    plotall(plot, source);

    // Use this to create output plots automatically
    //  d->SaveAsRoot("plotoptical.root");
    //  d->SaveAsPS("plotoptical.ps");

    return 1;
}

int plotoptical()
{
    return plotoptical("", "");
}