/*!\file PositiveDegreeDayx
 * \brief: create system matrices (stiffness matrix, loads vector)
 */

#include "./PositiveDegreeDayx.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../io/io.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

void PositiveDegreeDayx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){

// void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
//    note "v" prefix means 12 monthly means, ie time dimension
//    INPUT parameters: ni: working size of arrays
//    INPUT: surface elevation (m): hd(NA)
//    monthly mean surface sealevel temperature (degrees C): vTempsea(NA
//    ,NTIME) 
//    monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
//    OUTPUT: mass-balance (m/yr ice): agd(NA)
//    mean annual surface temperature (degrees C): Tsurf(NA)

  int    i, it, jj, itm;
  IssmDouble DT = 0.02, sigfac, snormfac;
  IssmDouble signorm = 5.5;      // signorm : sigma of the temperature distribution for a normal day 
  IssmDouble siglim;       // sigma limit for the integration which is equal to 2.5 sigmanorm
  IssmDouble signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
  IssmDouble siglimc, siglim0, siglim0c;
  IssmDouble tstep, tsint, tint, tstepc;
  int    NPDMAX = 1504, NPDCMAX = 1454;
  //IssmDouble pdds[NPDMAX]={0}; 
  //IssmDouble pds[NPDCMAX]={0};
  IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
  IssmDouble PDup, PDCUT = 2.0;    // PDcut: rain/snow cutoff temperature (C)
  IssmDouble tstar; // monthly mean surface temp
  
  IssmDouble* pdds=NULL; 
  IssmDouble* pds=NULL; 
  Element* element = NULL;
  
  pdds=xNew<IssmDouble>(NPDMAX+1); 
  pds=xNew<IssmDouble>(NPDCMAX+1); 
  
  // initialize PDD (creation of a lookup table)
  tstep = 0.1;
  tsint = tstep*0.5;
  sigfac = -1.0/(2.0*pow(signorm,2));
  snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
  siglim = 2.5*signorm;
  siglimc = 2.5*signormc;
  siglim0 =  siglim/DT + 0.5;
  siglim0c =  siglimc/DT + 0.5;
  PDup = siglimc+PDCUT;

  itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
  
  if (itm >= NPDMAX){
    _printLine_("increase NPDMAX in massBalance.cpp");
    exit (1);
      }
  for (it = 0; it < itm; it++){  
    //    tstar = REAL(it)*DT-siglim;
    tstar = it*DT-siglim;
    tint = tsint;
    pddt = 0.;
    for ( jj = 0; jj < 600; jj++){
      if (tint > (tstar+siglim)){break;}
      pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
      tint = tint+tstep;
    }
    pdds[it] = pddt*snormfac;
  }
  pdds[itm+1] = siglim + DT;
  
  //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
  tstepc = 0.1;
  tsint = PDCUT-tstepc*0.5;
  signormc = signorm - 0.5;
  sigfac = -1.0/(2.0*pow(signormc,2));
  snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
  siglimc = 2.5*signormc ;
  itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
  if (itm >= NPDCMAX){
    _printLine_("'increase NPDCMAX in p35com'");
    exit (1);
      }
  for (it = 0; it < itm; it++ ){
    tstar = it*DT-siglimc;
    //    tstar = REAL(it)*DT-siglimc;
    tint = tsint;          // start against upper bound
    pd = 0.;
    for (jj = 0; jj < 600; jj++){
      if (tint<(tstar-siglimc)) {break;}
      pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
      tint = tint-tstepc;
    }
    pds[it] = pd*snormfac;  // gaussian integral lookup table for snow fraction
  }
  pds[itm+1] = 0.;
  //     *******END initialize PDD
  
  for(i=0;i<elements->Size();i++){
    element=(Element*)elements->GetObjectByOffset(i);
    element->PositiveDegreeDay(pdds,pds,signorm);
  }
  /*free ressouces: */
  xDelete<IssmDouble>(pdds);
  xDelete<IssmDouble>(pds);
  
}
