/*
    Copyright 2012 Lumerical Solutions, Inc. All rights reserved.
*/
#include "magneticelectriclorentz.h"

using namespace std;
/*!
    \class MagneticElectricLorentzPlugin

    \brief This model implements both an electric and magnetic Lorentz material. It also
    includes constant electric and magnetic chi1
*/


const char* MagneticElectricLorentzPlugin::names[11] = {"\xcf\x87" "0 electric",              //0
                                                        "\xce\x94\xce\xb5 electric",        //1
                                                        "\xcf\x89" "0 electric",              //2
                                                        "\xce\xb4 electric",                //3
                                                        "\xcf\x87" "0 magnetic",              //4
                                                        "\xce\x94\xce\xbc magnetic",        //5
                                                        "\xcf\x89" "0 magnetic",              //6
                                                        "\xce\xb4 magnetic",                //7
                                                        "exclude electric",                 //8
                                                        "exclude magnetic",                 //9
                                                        0};

void MagneticElectricLorentzPlugin::initialize(const double** parameters, double dt)
{
    //loop over 3 axes
    for(int i=0; i<3; i++){

        chi1e[i] = float(parameters[0][i]);
        double deltaEps = parameters[1][i];
        double w0e = parameters[2][i];
        double de = parameters[3][i];
        chi1m[i] = float(parameters[4][i]);
        double deltaMu = parameters[5][i];
        double w0m = parameters[6][i];
        double dm = parameters[7][i];
        bool electricUpdate = !bool(parameters[8][i]);
        bool magneticUpdate = !bool(parameters[9][i]);

        //calculate update quantities
        if(electricUpdate) {
            e1[i] = float( (2.-w0e*w0e*dt*dt)/(de*dt+1.) );
            e2[i] = float( (de*dt-1.)/(de*dt+1.) );
            e3[i] = float( (deltaEps*w0e*w0e*dt*dt)/(de*dt+1.) );
        } else {
            e1[i] = e2[1] = e3[i] = chi1e[i] = 0.f;
        }
        if(magneticUpdate) {
            h1[i] = float( (2.-w0m*w0m*dt*dt)/(dm*dt+1.) );
            h2[i] = float( (dm*dt-1.)/(dm*dt+1.) );
            h3[i] = float( (deltaMu*w0m*w0m*dt*dt)/(dm*dt+1.) );
        } else {
            h1[i] = h2[i] = h3[i] = chi1m[i] = 0.f;
        }

    }
}


float MagneticElectricLorentzPlugin::calculateE(int i, float U, float V, float E,float *storage)
{
    float Pnm1 = storage[0];
    float Pn = storage[1];
    float Pnp1 = e1[i]*Pn + e2[i]*Pnm1 + e3[i]*E;

    storage[0] = Pn;
    storage[1] = Pnp1;

    return (V-Pnp1)/(U+chi1e[i]);

}

float MagneticElectricLorentzPlugin::calculateEx( float U, float V, float Ex, float* storage )
{
    return calculateE(0, U, V, Ex, storage);
}

float MagneticElectricLorentzPlugin::calculateEy( float U, float V, float Ey, float* storage )
{
    return calculateE(1, U, V, Ey, storage);
}

float MagneticElectricLorentzPlugin::calculateEz( float U, float V, float Ez, float* storage )
{
    return calculateE(2, U, V, Ez, storage);
}

float MagneticElectricLorentzPlugin::calculateH( int i, float U, float V, float H, float *storage )
{
    float Mnm1 = storage[0];
    float Mn = storage[1];
    float Mnp1 = h1[i]*Mn + h2[i]*Mnm1 + h3[i]*H;

    storage[0] = Mn;
    storage[1] = Mnp1;

    return (V-Mnp1)/(U+chi1m[i]);
}

float MagneticElectricLorentzPlugin::calculateHx( float U, float V, float Hx, float* storage )
{
    return calculateH(0, U, V, Hx, storage);
}

float MagneticElectricLorentzPlugin::calculateHy( float U, float V, float Hy, float* storage )
{
    return calculateH(1, U, V, Hy, storage);
}

float MagneticElectricLorentzPlugin::calculateHz( float U, float V, float Hz, float* storage )
{
    return calculateH(2, U, V, Hz, storage);
}

MATERIAL_PLUGIN(MagneticElectricLorentzPlugin);