solar_panel.cc

#include "solar_panel.h"
#include <string.h>

// Prototypes
ml_matrix sun_vector(double jd);
ml_matrix earth_rot(double jcent);
double gms_time(double jd);
double jd_to_midnight(double jd);
double jd_to_jcent(double jd);

// The number of radians in 1 degree.
const double DEGREES_TO_RADIANS = 0.0174532925199433;
// The number of degrees in 1 revolution.
const double REVS_TO_DEGREES = 360.0;
// Julian date value for 2000-01-01 12:00:00+00.
const double JD_2000 = 2451545.0;
// Inverse of DAYS_TO_SECS.
const double SECS_TO_DAYS = 1.157407407407407e-05;
//å Inverse of CENTURIES_TO_DAYS.
const double DAYS_TO_CENTURIES = 2.737850787132101e-05;



solar_panel :: solar_panel(dsim_model_setup *setup) : dsim_model(setup)
{
}

void solar_panel :: initialize_data()
{
        dsim_model::initialize_data();
        
        double power_total = 0;
        ml_matrix normal(3,1);          normal(1,1) = 1;
        ml_matrix area(1,1);            area(1,1) = 25;
        ml_matrix efficiency(1,1);      efficiency(1,1) = 0.22;
    
        power_dsim                      =       create_output("power",                  sd_type_double, &power_total,   1,1,    "W",    "Power panel"); 
        normal_dsim                     =       create_parameter("normal",      sd_type_matrix, &normal,                3,1,    "",             "Solar panel cell face normal");        
        area_dsim                       =       create_parameter("area",        sd_type_matrix, &area,                  1,1,    "m^2",  "Solar panel area");    ;
        efficiency_dsim         =       create_parameter("efficiency",  sd_type_matrix, &efficiency,    1,1,    "",             "Solar panel efficiency");      
        
        power_dsim.mark_telemetry();
        
        configure_timestep(true, false, false);
        
}

void solar_panel :: initialization_complete()
{
    dsim_model::initialization_complete();
}

void solar_panel :: pre_calculate()
{
        dsim_model::initialize_timestep();      
        
        double jd               = sim_jd();
    ml_matrix solar_flux    = 1000*earth_rot(jd)*sun_vector(jd);
        ml_matrix n             = normal_dsim.value_as_matrix();
        ml_matrix area          = area_dsim.value_as_matrix();
        ml_matrix eff           = efficiency_dsim.value_as_matrix();

        double p;
        double power_total = 0;

        for( int k = 1; k <= n.cols(); k++ )
        {
                p = area.get(1,k)*eff.get(1,k)*solar_flux.dot(n.get_col(k));
                if( p < 0 ) p = 0;
                power_total += p;
        }
        
        power_dsim.set_value_as_double( power_total );
}

double jd_to_midnight(double jd)
{
    double jd_noon = floor(jd);
    return (jd - jd_noon >= 0.5) ? (jd_noon + 0.5) : (jd_noon - 0.5);
}

double jd_to_jcent(double jd)
{
        return (jd - JD_2000) * DAYS_TO_CENTURIES;
}

double gms_time(double jd)
{
    // Julian days at 0h UT
    double jd_mid = jd_to_midnight(jd);
    double t_diff = jd_to_jcent(jd_mid);
    
    // This organization maximizes the precision
    double gmst = ((0.093104 - 6.2e-6*t_diff)*t_diff + 8640184.812866)*t_diff + 24110.54841;
    
    // Account for earth rotation
    gmst *= SECS_TO_DAYS;
    gmst += jd-jd_mid;
    
    // Limit to the range of 0 to 360 degrees
    gmst = fmod(gmst*360, 360); 
    if (gmst < 0)
    {
        gmst = gmst + 360;
    }
    return gmst; 
}  


ml_matrix earth_rot(double jd)
{
    // Find Greenwich Mean Sidereal Time
    double gmst = gms_time(jd);
    double sin_gmst = sin(gmst * DEGREES_TO_RADIANS);
    double cos_gmst = cos(gmst * DEGREES_TO_RADIANS);
    ml_matrix b(3,3);
    b(1,1) = cos_gmst;
    b(1,2) = sin_gmst;
    b(2,1) = -sin_gmst;
    b(2,2) = cos_gmst;
    b(3,3) = 1;
    return b;
}


ml_matrix sun_vector(double jd)
{
    // Mean anomaly
    double t_diff = jd - JD_2000;
    double g = fmod((357.528 + 0.9856003 * t_diff),REVS_TO_DEGREES) * DEGREES_TO_RADIANS; 
    // Ecliptic longitude
    double lam = (fmod(280.460 + 0.9856474 * t_diff, REVS_TO_DEGREES) + 1.915*sin(g) + 0.02*sin(2.0*g)) * DEGREES_TO_RADIANS;
    
    // Obliquity of the ecliptic
    double ob_of_e = (23.439 - 4.00e-7 * t_diff) * DEGREES_TO_RADIANS;
    
    // Equatorial rectangular coordinates of the Sun
    double s_lam  = sin(lam);
    ml_matrix s_vec(3,1);
    s_vec(1,1) = cos(lam);
    s_vec(2,1) = cos(ob_of_e)*s_lam;
    s_vec(3,1) = sin(ob_of_e)*s_lam;
    return s_vec;
}



extern "C"
{
    
    dsim_model *SolarPanelModel(void *setup);

        dsim_model *SolarPanelModel(void *setup)
        {
                return new solar_panel((dsim_model_setup *)setup);
        }
}