/*Project leader: Marco Montecchi
             ENEA (Italy)
             email: marco.montecchi@enea.it

Algorithm3 to calculate the sun position described in Solar Energy 86 (2012) 1323-1337
            Roberto Grena
            ENEA (Italy)
            email: roberto.grena@enea.it

Italian Solar Radiation Atlas http://www.solaritaly.enea.it/indexEn.php
            Irena Balog and Giampaolo Caputo
            ENEA (Italy)
            email: irena.balog@enea.it  giampaolo.caputo@enea.it

Porting to Windows
             Alberto Mittiga
             ENEA (Italy)
             email: alberto.mittiga@enea.it


MeanPTC is a calculator of incidence and acceptance angles mean values, weighthed on DNI data.
MeanPTC includes an alghoritm for calculating the Sun position given the site coordinates and time,
    thus only DNI data are mandatory

   Copyright (C) 2023  Marco Montecchi

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/
#include "meanptc.h"
#include "ui_meanptc.h"
#include <iostream>
#include <math.h>
#include <QFileDialog>
#include <QtWidgets>
#include <qwt_plot.h>
#include <qwt_plot_curve.h>
#include <qwt_series_data.h>
#include <qwt_plot_grid.h>
#include <QMessageBox>

using namespace std;

//global variables
int nLFCrows;
double deg2rad=acos(-1.)/180.;
QString curDir,fileDNI[10];
double focal,radius;
double Rheight,Raperture,PMwidth;
double ThetaL,ThetaT,ThetaLS,Xmin,Xmax,Orientamento;
double meanTheta,stdTheta,meanPhi,stdPhi;
double PI=acos(-1);
double PI2=2.*PI;
double PIM=0.5*PI;
double ThetaDist[90];//distribution of Theta
double PhiDist[100];//distribution of Phi
double LatLongTdel[10][3];

// input data:
double UT;
int Day;
int Month;
int Year;
double Dt;
double Longitude;
double Latitude;
double Pressure;
double Temperature;

//output data
double RightAscension;
double Declination;
double HourAngle;
double Zenith;
double Azimuth;
double uS[3];

//functions



MeanPTC::MeanPTC(QWidget *parent)
    : QMainWindow(parent)
    , ui(new Ui::MeanPTC)
{
    ui->setupUi(this);

    connect(ui->pushButton_SunPos,SIGNAL(clicked()),this, SLOT(computingSunPos())   );
    connect(ui->pushButton_compute,SIGNAL(clicked()),this, SLOT(computing())        );
    connect(ui->pushButton_load,SIGNAL(clicked()),this, SLOT(computingYear())       );
    connect(ui->pushButton_contributors,SIGNAL(clicked()),this, SLOT(contributors()));
    connect(ui->comboBox_site,SIGNAL(currentIndexChanged(int)),this,SLOT(setZenAzi()));

    // parameter initialization
#ifdef __unix__
#define IS_POSIX 1
#else
#define IS_POSIX 0
#endif

    QDir dir;  //current directory
    dir.cdUp();//cd ..
    dir.cdUp();//cd ..
    if (IS_POSIX == 1) {
        //Linux path initialization
        //nothing to do
    }
    else {
        //windows path inizialization
        dir.cdUp();//cd ..
    }
    curDir=dir.absolutePath();
    QString fSiteCoord=curDir+"/qtSource/MeanPTC/SiteCoord.txt";
    curDir=curDir+"/DNIdata";
    printf("Curdir= %s \n",curDir.toStdString().c_str());
    //building site coordinates combobox
    QFile file0(fSiteCoord);
    if (!file0.open (QIODevice::ReadOnly | QIODevice::Text))
        return;
    QTextStream stream ( &file0 );
    QString line;
    ui -> comboBox_site-> clear();
    ui -> comboBox_site->addItem("Custom");
    int iSite=0;
    do{
        line = stream.readLine();
        QString label=line.simplified();
        ui -> comboBox_site->addItem(label);
        line=stream.readLine();
        line=line.simplified();
        QStringList List;
        List =line.split(" ");
        int nV=List.count();
        //printf("nV=%d\n",nV);
        if(nV!=3)
            return;
        QString pezzo;
        pezzo=List.at(0).toLocal8Bit().constData();
        LatLongTdel[iSite][0]=pezzo.toDouble();
        pezzo=List.at(1).toLocal8Bit().constData();
        LatLongTdel[iSite][1]=pezzo.toDouble();
        pezzo=List.at(2).toLocal8Bit().constData();
        LatLongTdel[iSite][2]=pezzo.toDouble();
        line=stream.readLine();
        fileDNI[iSite]=line;
        printf("Site =%s Latitude=%f Longitude=%f with fileDNI=%s and Tdelay=%f is added\n",
               label.toStdString().c_str(),LatLongTdel[iSite][0],LatLongTdel[iSite][1],line.toStdString().c_str(),LatLongTdel[iSite][2]);
        iSite++;
    }while(!stream.atEnd());
    file0.close();
    fflush(stdout);
    setZenAzi();
}

MeanPTC::~MeanPTC()
{
    delete ui;
}

void MeanPTC::contributors(){
    QMessageBox msgBox;
    QString msg="Contributors:\nMarco Montecchi (project leader)\nRoberto Grena (Sun position algorithm)\nIrena Balog and Giampaolo Caputo (Italian Solar Radiation Atlas)\nAlberto Mittiga (porting to Windows)\n";
    msg=msg+"\nLicense:\nThe code source files are distributed as open source software under the GNU General Public License as published by the Free Software Foundation version 3";
    msgBox.setText(msg);
    msgBox.setStandardButtons(QMessageBox::Ok);
    msgBox.exec();
    return;
}


void MeanPTC::setZenAzi(){
    int iSite=ui -> comboBox_site->currentIndex();
    if(iSite==0){
        ui->groupBox_3->setEnabled(true);
    }
    else{
        ui->doubleSpinBox_Lat->setValue(LatLongTdel[iSite-1][0]);
        ui->doubleSpinBox_Long->setValue(LatLongTdel[iSite-1][1]);
        ui->doubleSpinBox_Tdelay->setValue(LatLongTdel[iSite-1][2]);
        ui->groupBox_3->setEnabled(false);
    }
}


void MeanPTC::computingSunPos(){
    printf("-> computing Sun position at ");
    Orientamento=ui->doubleSpinBox_Orientamento->value();
    Orientamento=Orientamento*deg2rad;
    Latitude=ui->doubleSpinBox_Lat->value();
    Longitude=ui->doubleSpinBox_Long->value();
    Latitude=Latitude*deg2rad;
    Longitude=Longitude*deg2rad;
    QDate date=ui->dateTimeEdit->date();
    Year=date.year();
    Month=date.month();
    Day=date.day();
    QTime time=ui->dateTimeEdit->time();
    int hou=time.hour();
    int min=time.minute();
    int sec=time.second();
    UT=double(hou)+double(min)/60.+double(sec)/3600.;
    printf("Year=%d Month=%d Day=%d Hour=%f\n",Year,Month,Day,UT);
    Algorithm3();
    thetasANDuS();
}


void MeanPTC::computing(){
    printf("-> computing at a given time\n");
    //reset counters
    for(int i=0;i<90;i++)
        ThetaDist[i]=0.;
    for(int i=0;i<100;i++)
        PhiDist[i]=0.;
    computingSunPos();
    if(Zenith/deg2rad>90.){
        printf("Zenith>90° -> abort\n");
        return;
    }

    int i0PT2LF=ui->tabWidget->currentIndex();
    if(i0PT2LF==0){//PTC
        focal=ui->doubleSpinBox_focal->value();
        radius=ui->doubleSpinBox_radius->value();
        Xmin=ui->doubleSpinBox_Xmin->value();
        Xmax=ui->doubleSpinBox_Xmax->value();
        kernelComputing(1.);
    }
    else{//LFC
        nLFCrows=ui->spinBox_Nrows->value();
        Rheight=ui->doubleSpinBox_LFCheight->value();
        Raperture=ui->doubleSpinBox_LFCaperture->value();
        PMwidth=ui->doubleSpinBox_LFCwidth->value();
        kernelComputing2(1.);
    }
    Averagetor();
    ui->lineEdit_Theta->setText(QString::number(meanTheta,'f',1)+" +- "+QString::number(stdTheta,'f',1));
    ui->lineEdit_Phi->setText(QString::number(meanPhi,'f',1)+" +- "+QString::number(stdPhi,'f',1));
}



void MeanPTC::computingYear(){
    //reset counters
    for(int i=0;i<90;i++)
        ThetaDist[i]=0.;
    for(int i=0;i<100;i++)
        PhiDist[i]=0.;

    int Ncol,NcolDate,NcolEle,NcolAzi,NcolDNI,iSep;
    double weight,Tdelay;
    int iSite=ui->comboBox_site->currentIndex();
    int i0PT2LF=ui->tabWidget->currentIndex();
    int iE0Z1=ui->comboBox_E0Z1->currentIndex();
    Latitude=ui->doubleSpinBox_Lat->value();
    Longitude=ui->doubleSpinBox_Long->value();
    Latitude=Latitude*deg2rad;
    Longitude=Longitude*deg2rad;
    Orientamento=ui->doubleSpinBox_Orientamento->value();
    Orientamento=Orientamento*deg2rad;
    if(i0PT2LF==0){//PTC
        focal=ui->doubleSpinBox_focal->value();
        radius=ui->doubleSpinBox_radius->value();
        Xmin=ui->doubleSpinBox_Xmin->value();
        Xmax=ui->doubleSpinBox_Xmax->value();
    }
    else{//LFC
        nLFCrows=ui->spinBox_Nrows->value();
        Rheight=ui->doubleSpinBox_LFCheight->value();
        Raperture=ui->doubleSpinBox_LFCaperture->value();
        PMwidth=ui->doubleSpinBox_LFCwidth->value();
    }
    if(iSite==0){
        Ncol=ui->spinBox_Ncol->value();
        NcolDate=ui->spinBox_ColDate->value();
        NcolEle=ui->spinBox_ColEle->value();
        NcolAzi=ui->spinBox_ColAzi->value();
        NcolDNI=ui->spinBox_ColDNI->value();
        iSep=ui->comboBox_SepChar->currentIndex();
    }
    else{
        Ncol=3;
        NcolDate=1;
        NcolDNI=3;
        iSep=0;
    }
    Tdelay=ui->doubleSpinBox_Tdelay->value();
    QString sep;
    if(iSep==0)
        sep=" ";
    else
        sep=";";
    QString line;
    QString fileName;
    if(iSite!=0)
        fileName=curDir+"/"+fileDNI[iSite-1];
    else
        fileName= QFileDialog::getOpenFileName(
        this,
        "Choose a Yearly Data File", //window title
        curDir,                 //initial directory
        " (*.*)"); //file extension
    //printf("fileName= %s\n",fileName.toStdString().c_str());
    QFile file(fileName);
    if (!file.open (QIODevice::ReadOnly | QIODevice::Text))
        return;
    QTextStream stream ( &file );
    do{
        line=stream.readLine();
        line=line.simplified();
        //printf("%s\n",line.toStdString().c_str());
        QStringList List;
        List =line.split(sep);
        int nV=List.count();
        //printf("nV=%d\n",nV);
        if(nV!=Ncol)
            continue;
        QString pezzo;
        pezzo=List.at(NcolDNI-1).toLocal8Bit().constData();
        weight=pezzo.toDouble()/1000.;
        //printf("weight= %f\n",weight);
        if(weight>0.){
            QDate date = QDate::fromString(List.at(NcolDate-1),"yyyy-MM-dd");
            Year=date.year();
            Month=date.month();
            Day=date.day();
            QTime time=QTime::fromString(List.at(NcolDate),"HH:mm");
            int hou=time.hour();
            int min=time.minute();
            int sec=time.second();
            UT=double(hou)+double(min)/60.+double(sec)/3600.-Tdelay/60.;
            printf("Year=%d Month=%d Day=%d Hour=%f weight=%f\n",Year,Month,Day,UT,weight);
            if(iSite!=0 ||NcolEle==0 || NcolAzi==0){
                Algorithm3();
            }
            else{
                pezzo=List.at(NcolAzi-1).toLocal8Bit().constData();
                Azimuth=pezzo.toDouble();
                pezzo=List.at(NcolEle-1).toLocal8Bit().constData();
                if(iE0Z1==0)
                    Zenith=90.-pezzo.toDouble();
                else
                    Zenith=pezzo.toDouble();
//                printf("Zenith=%f Azimuth=%f\n",Zenith,Azimuth);
                Azimuth=Azimuth*deg2rad;
                Zenith=Zenith*deg2rad;
            }
            thetasANDuS();
            if(Zenith>90./deg2rad)
                continue;
            if(i0PT2LF==0){//PTC
                weight=weight*cos(ThetaLS);
                kernelComputing(weight);
            }
            else//LF
                kernelComputing2(weight);
        }
        fflush(stdout);
    }while(!stream.atEnd());
    file.close();
    Averagetor();
    ui->lineEdit_Theta->setText(QString::number(meanTheta,'f',1)+" +- "+QString::number(stdTheta,'f',1));
    ui->lineEdit_Phi->setText(QString::number(meanPhi,'f',1)+" +- "+QString::number(stdPhi,'f',1));
}


void MeanPTC::kernelComputing(double weight){
    //Parabolic-Trough Collector
    int Nstep=101;
    double Xstep=(Xmax-Xmin)/(Nstep-1);
    double x,y,alpha,costheta,theta,t,phi;//,uSpt[3],uN[3],uR[3]
    for(int i=0;i<Nstep;i++){
        x=Xmin+Xstep*i;
        y=0.25/focal*x*x;
        alpha=atan(0.5/focal*x);
        costheta=cos(ThetaLS)*cos(alpha);
        theta=acos(costheta);
        t=(focal-y)/(2.*costheta*cos(alpha)-cos(ThetaLS));
//        uSpt[0]=0.;
//        uSpt[1]=sin(ThetaLS);
//        uSpt[2]=cos(ThetaLS);
//        uN[0]=-sin(alpha);
//        uN[1]=0.;
//        uN[2]=cos(alpha);
//        for(int k=0;k<3;k++)
//            uR[k]=2.*costheta*uN[k]-uSpt[k];
//        printf("uR: %f, %f, %f\n",uR[0],uR[1],uR[2]);
//        printf("t*uR[2]+y: %f\n",t*uR[2]+y);
        phi=radius/t*1000.;//mrad
        //printf("x=%f Theta= %f t=%f Phi=%f\n",x,theta/deg2rad,t,phi);
//        fflush(stdout);
        int iT=int(theta/deg2rad);
        if(iT<90)
            ThetaDist[iT]=ThetaDist[iT]+weight;
        int iP=int(phi);
        if(iP<100)
            PhiDist[iP]=PhiDist[iP]+weight;
    }
}

void MeanPTC::kernelComputing2(double weight){
    //Linear Fresnel Collector
    double MirWidth=PMwidth/double(nLFCrows);//width of the single mirror row
    double Xstep=(PMwidth-MirWidth)/double(nLFCrows-1);//step between joint-pivots
    double x,alpha,beta,theta,t,phi,uN[3];//,uR[3]
    for(int i=0;i<nLFCrows;i++){
        x=-PMwidth/2.+MirWidth/2.+i*Xstep;//joint-pivot position of mirror i-row
        alpha=atan(x/Rheight);//angle of the straight line crossing x and the aperture centre
        beta=ThetaT-(ThetaT-alpha)/2.;//tilt angle of mirror i-row
        uN[0]=-sin(beta);
        uN[1]=0.;
        uN[2]=cos(beta);
        //incidence angle from scalar product Sun-ray and mirror-normal
        double prosca=0;
        for(int k=0;k<3;k++)
            prosca=prosca+uN[k]*uS[k];
        theta=acos(prosca);
        //for(int k=0;k<3;k++)
        //    uR[k]=2.*prosca*uN[k]-uS[k];
        t=sqrt(x*x+Rheight*Rheight)/cos(ThetaLS);//path length to the receiver
        //printf("t*uR: %f, %f, %f\n",t*uR[0],t*uR[1],t*uR[2]);
        phi=0.5*Raperture*cos(alpha)/t*1000.;//mrad
        //printf("x=%f alpha=%f beta=%f theta=%f t=%f phi=%f\n",x,alpha/deg2rad,beta/deg2rad,theta/deg2rad,t,phi);
        int iT=int(theta/deg2rad);
        if(iT<90)
            ThetaDist[iT]=ThetaDist[iT]+weight*cos(theta);
        int iP=int(phi);
        if(iP<100)
            PhiDist[iP]=PhiDist[iP]+weight*cos(theta);
    }
    //fflush(stdout);
}

void MeanPTC::Averagetor(){
    double SumT=0.,SumTT=0.,SumP=0.,SumPP=0.,SumW=0.;
    for(int i=1;i<90;i++){
        SumW=SumW+ThetaDist[i];
        SumT=SumT+double(i)*ThetaDist[i];
        SumTT=SumTT+double(i*i)*ThetaDist[i];
    }
    for(int i=1;i<100;i++){
        SumP=SumP+double(i)*PhiDist[i];
        SumPP=SumPP+double(i*i)*PhiDist[i];
    }
    meanTheta=SumT/SumW;
    stdTheta=sqrt(SumTT/SumW-meanTheta*meanTheta);
    meanPhi=SumP/SumW;
    stdPhi=sqrt(SumPP/SumW-meanPhi*meanPhi);
    printf("->Averagetor: \n\tTheta= %f +- %f \n\tPhi= %f +- %f\n",meanTheta,stdTheta,meanPhi,stdPhi);
    //saving distributions
    QString fileDistT=curDir+"/ThetaDistribution.txt";
    QString fileDistP=curDir+"/PhiDistribution.txt";
//    return;
    QFile fileDT(fileDistT);
    if (!fileDT.open(QIODevice::WriteOnly | QIODevice::Text)){
        printf("error on opening file =%s \n", fileDistT.toStdString().c_str());
        fflush(stdout);
        return;
    }
    QTextStream out(&fileDT);
    out << "Theta(deg)" << "\t"<<"frequency\n";
    for(int i=0;i<90;i++)
        out << i << "\t"<<ThetaDist[i]/SumW<<"\n";
    fileDT.close();
    QFile fileDP(fileDistP);
    if (!fileDP.open(QIODevice::WriteOnly | QIODevice::Text)){
        printf("error on opening file =%s \n", fileDistP.toStdString().c_str());
        fflush(stdout);
        return;
            }
    QTextStream out2(&fileDP);
    out2 << "Phi(mrad)" << "\t"<<"frequency\n";
    for(int i=0;i<100;i++)
        out2 << double(i) << "\t"<<PhiDist[i]/SumW<<"\n";
    fileDP.close();

    //theta graphs
    double xPlot[100],yPlot[100];
    for(int i=0;i<90;i++){
        xPlot[i]=double(i);
        yPlot[i]=ThetaDist[i]/SumW;
    }
    QwtPlot *Gt=new QwtPlot();
    QwtPlotCurve *curve1=new QwtPlotCurve("Curve 1");
    curve1 -> setSamples(xPlot,yPlot,90);
    curve1 -> setPen(QPen(Qt::blue,2,Qt::SolidLine));
    curve1 -> attach(Gt);
    Gt -> setAxisTitle(0,"Frequency");
    Gt -> setAxisTitle(2,"Theta (deg)");
    Gt -> setAxisScale(2,0,90,0);
    Gt -> setAutoReplot();
    Gt -> show();

    //phi graph
    double Phimax=0.;
    for(int i=0;i<100;i++){
        xPlot[i]=double(i);
        yPlot[i]=PhiDist[i]/SumW;
        if(yPlot[i]>0.)
            Phimax=double(i);
    }
    Phimax++;
    QwtPlot *Gp=new QwtPlot();
    QwtPlotCurve *curve2=new QwtPlotCurve("Curve 2");
    curve2 -> setSamples(xPlot,yPlot,100);
    curve2 -> setPen(QPen(Qt::red,2,Qt::SolidLine));
    curve2 -> attach(Gp);
    Gp -> setAxisTitle(0,"Frequency");
    Gp -> setAxisTitle(2,"Phi (mrad)");
    Gp -> setAxisScale(2,0,Phimax,0);
    Gp -> setAutoReplot();
    Gp -> show();
}


void MeanPTC::Algorithm3() {
    //auxiliary
    double t, te, wte, cl,sl,ce,se,epsi,
        sp, cp, sd, cd, sH, cH, se0, ep, De,lambda;
    int yt, mt;

    if (Month <= 2) {
        mt = Month + 12;
        yt = Year - 1;
    } else {
        mt = Month;
        yt = Year;
    }

    // Montecchi's addictions
    Dt=96.4+0.567*(yt-2061);
    Pressure=1;//atmosphere
    Temperature=20.;//°C

    t = double(int(365.25*double(yt-2000)) + int(30.6001*double(mt+1)) - int(0.01*double(yt)) + Day) + 0.0416667*UT - 21958.0;
    te = t + 1.1574e-5*Dt;

    wte = 0.0172019715*te;

    lambda = -1.388803 + 1.720279216e-2*te + 3.3366e-2*sin(wte - 0.06172) + 3.53e-4*sin(2.0*wte - 0.1163);

    epsi = 4.089567e-1 - 6.19e-9*te;

    sl = sin(lambda);
    cl = cos(lambda);
    se = sin(epsi);
    ce = sqrt(1-se*se);

    RightAscension = atan2(sl*ce, cl);
    if (RightAscension < 0.0)
        RightAscension += PI2;

    Declination = asin(sl*se);

    HourAngle = 1.7528311 + 6.300388099*t + Longitude - RightAscension;
    HourAngle = fmod(HourAngle + PI, PI2) - PI;
    if (HourAngle < -PI) HourAngle += PI2;

    //if (flag == 's') return;

    sp = sin(Latitude);
    cp = sqrt((1-sp*sp));
    sd = sin(Declination);
    cd = sqrt(1-sd*sd);
    sH = sin(HourAngle);
    cH = cos(HourAngle);
    se0 = sp*sd + cp*cd*cH;
    ep = asin(se0) - 4.26e-5*sqrt(1.0-se0*se0);
    Azimuth = atan2(sH, cH*sp - sd*cp/cd);

    if (ep > 0.0)
        De = (0.08422*Pressure) / ((273.0+Temperature)*tan(ep + 0.003138/(ep + 0.08919)));
    else
        De = 0.0;

    Zenith = PIM - ep - De;
//    printf("Algorithm3: t=%f Zenith=%f Azimuth=%f\n",t,Zenith/deg2rad,Azimuth/deg2rad);
}

void  MeanPTC::thetasANDuS(){
    ThetaL=atan(cos(Azimuth-Orientamento)*tan(Zenith));
    ThetaT=atan(sin(Azimuth-Orientamento)*tan(Zenith));
    ThetaLS=asin(cos(Azimuth-Orientamento)*sin(Zenith));
    uS[0]=-sin(Azimuth-Orientamento)*sin(Zenith);
    uS[1]=-cos(Azimuth-Orientamento)*sin(Zenith);
    uS[2]=cos(Zenith);
    ui->lineEdit_Azi->setText(QString::number(Azimuth/deg2rad,'f',1));
    ui->lineEdit_Elev->setText(QString::number(Zenith/deg2rad,'f',1));
    ui->lineEdit_ThetaL->setText(QString::number(ThetaL/deg2rad,'f',1));
    ui->lineEdit_ThetaT->setText(QString::number(ThetaT/deg2rad,'f',1));
    ui->lineEdit_ThetaLS->setText(QString::number(ThetaLS/deg2rad,'f',1));
//    printf("Zenith=%f Azimuth=%f ThetaL=%f ThetaT=%f ThetaLS=%f uS=(%f, %f, %f)\n",
//           Zenith/deg2rad,Azimuth/deg2rad,ThetaL/deg2rad,ThetaT/deg2rad,ThetaLS/deg2rad,uS[0],uS[1],uS[2]);
}
