/**

    @class Assemble_ZPT
    @author Brian Magill
    @creationdate 8/23/2006

    $Date:$
    $Revision:$

    @copyright
     (©) Copyright 2006 by GATS Inc.
     11864 Canon Blvd., Suite 101, Newport News, VA 23606

     All Rights Reserved. No part of this software or publication may be
     reproduced, stored in a retrieval system, or transmitted, in any form
     or by any means, electronic, mechanical, photocopying, recording, or
     otherwise without the prior written permission of GATS Inc.
    
    @brief Assembles ZPT profile from time, Ephemeris, Solar Flux, and NCEP data  
    
*/

#include <iostream>

#include "Assemble_ZPT.h"
#include "ZT_Profile.h"
#include "PressProfile.h"
#include "MSISCalc.h"
#include "SofieZTMerge.h"
#include "Create_ZPT.h"
#include "MiscExceptions.h"
#include "Vec2Val.hpp"
#include "MSISGases.h"
#include "MSISXRatio.h"

using namespace std;

Assemble_ZPT::Assemble_ZPT(long date, double sec, SolarFlux const & flux, EphemData const & ephemstd,
                 valarray<double> const & z_ncep, valarray<double> const & p_ncep, 
                 valarray<double> const & t_ncep)
    {
    
    const unsigned long N = 76;
    const double alt_grid = 2.0;
    unsigned long i;
    double alt;
    double tem;
    
    vector<double> Zout;     
    vector<double> Tout; 
    vector<double> Pout; 
    
    vector<double> z_msis;
    vector<double> t_msis;
    vector<double> O_msis;   //  atomic oxygen's mixing ratio
    vector<double> O2_msis;  //  diatomic oxygen's mixing ratio
    vector<double> N2_msis;  //  diatomic nitrogen's mixing ratio

//  Before going any further, make sure that the NCEP data is "reasonable"         
//    Check_NCEP(z_ncep, p_ncep, t_ncep);

    PressProfile ncepPressProf(z_ncep, p_ncep) ;
    
    ZT_Profile ncepZTProf(z_ncep, t_ncep);
    ZT_Profile msisZTProf;
    ZT_Profile mergedZTProf;

    MSISCalc msisCalc(date, flux, ephemstd);
    MSISGases msisDensities(date, flux, ephemstd);
    MSISXRatio xRatio;
        
    SofieZTMerge mergeProfiles;
        
//  Set up MSIS temperature profile     
//  Also calculate MSIS atomic O, O2, and N2 profiles

    for(i = 0; i < N; i++)
        {
        alt = alt_grid*i;
        z_msis.push_back(alt);
        tem = msisCalc(static_cast<float>(alt), 
                        static_cast<float>(sec));
        t_msis.push_back(tem);
        
        xRatio = MSISXRatio(msisDensities(static_cast<float>(alt), static_cast<float>(sec)) );
        O_msis.push_back(static_cast<double>(xRatio.getO_xratio() ) );
        O2_msis.push_back(static_cast<double>(xRatio.getO2_xratio() ) );
        N2_msis.push_back(static_cast<double>(xRatio.getN2_xratio() ) );
        
        }

    
    msisZTProf = ZT_Profile(z_msis, t_msis);
    
    mergedZTProf = mergeProfiles(ncepZTProf, msisZTProf);
    
    Create_ZPT createZPT(mergedZTProf, ncepPressProf, ephemstd);

//  Convert output to valarrays

    Vec2Val (createZPT.getAltitude(), altitude);
    Vec2Val (createZPT.getTemperature(), temperature);
    Vec2Val (createZPT.getPressure(), pressure);

    Vec2Val (z_msis,  xRatioAlt);
    Vec2Val (O_msis,  xRatio_O); 
    Vec2Val (O2_msis, xRatio_O2); 
    Vec2Val (N2_msis, xRatio_N2); 
    
    };                 

void Assemble_ZPT::getAltitude(std::valarray<double>  & z ) 
    { 
    z.resize(altitude.size() );
    z = altitude; 
    }; 

void Assemble_ZPT::getTemperature(std::valarray<double>  & tem ) 
    { 
    tem.resize(temperature.size() );
    tem = temperature;
    }; 

void Assemble_ZPT::getPressure(std::valarray<double>  & press) 
    { 
    press.resize(pressure.size() );
    press = pressure; 
    }; 

void Assemble_ZPT::getXRatioAlt(std::valarray<double>  & z )
    { 
    z.resize(xRatioAlt.size() );
    z = xRatioAlt; 
    }; 

void Assemble_ZPT::getxRatio_O(std::valarray<double>  & specie ) 
    { 
    specie.resize(xRatio_O.size() );
    specie = xRatio_O; 
    }; 

void Assemble_ZPT::getxRatio_O2(std::valarray<double>  & specie ) 
    { 
    specie.resize(xRatio_O2.size() );
    specie = xRatio_O2; 
    }; 

void Assemble_ZPT::getxRatio_N2(std::valarray<double>  & specie ) 
    { 
    specie.resize(xRatio_N2.size() );
    specie = xRatio_N2; 
    }; 

/**
  @def Check_NCEP(valarray<double> const & z_ncep, 
                     valarray<double> const & p_ncep, 
                     valarray<double> const & t_ncep)
  
  
 
  @param[in] z_ncep  NCEP altitude array
  @param[in] p_ncep  NCEP pressure array  
  @param[in] t_ncep  NCEP temperature array  
 
  @brief Checks incoming NCEP arrays for errors that might crash the program
 
 */

void Assemble_ZPT::Check_NCEP(valarray<double> const & z_ncep, 
                     valarray<double> const & p_ncep, 
                     valarray<double> const & t_ncep)
    {
    
    string errMsg;
/*
for(int i = 0 ; i< (int)z_ncep.size() ; ++i ) 
	std::cout << z_ncep[i] << "  " << p_ncep[i] << "  " << t_ncep[i] << std::endl;
exit(23);
*/

    if (z_ncep.size() < 2 )
        {
        errMsg = "Assemble_ZPT: NCEP arrays do not have enough data";
        THROW_GATS_EXCEPTION(BadDataException, errMsg);
        }
    
    if (z_ncep.size() != p_ncep.size() )
        {
        errMsg = "Assemble_ZPT: NCEP altitude array's size does not match pressure array's";
        THROW_GATS_EXCEPTION(BoundaryException, errMsg);        
        }
        
    if (z_ncep.size() != t_ncep.size() )
        {
        errMsg = "Assemble_ZPT: NCEP altitude array's size does not match pressure array's";   
        THROW_GATS_EXCEPTION(BoundaryException, errMsg);     
        }
    
    for(unsigned long i = 0; i < z_ncep.size(); i++)
        {
        
        if (p_ncep[i] <= 0.0)
            {
            errMsg = "Assemble_ZPT: NCEP pressure array contains nonphysical values!";
            THROW_GATS_EXCEPTION(BadDataException, errMsg);
            } 
            
        if (t_ncep[i] <= 0.0)
            {
            errMsg = "Assemble_ZPT: NCEP temperature array contains nonphysical values!";
            THROW_GATS_EXCEPTION(BadDataException, errMsg);
            }                       
        }
     
    };