// $Id: scene.h,v 1.5 2010/09/17 15:22:56 thernis Exp $

#ifndef SCENE_H
#define SCENE_H

#include <string>
#include "camera.h"
#include "sun.h"
#include "Cbasis.h"
#include "Cvec.h"
#include "Clos.h"
#include "CModelBase.h"
#include "rtmiscfunc.h"
#include "ModelPosition.h"
#include "physicsbase.h"

//! Defines the Scene: Camera + Sun + Density
class Scene
{
  
  protected:
    Cbasis abs,obs; //!< absolute and observer coordinate system
    CModelBase *pmod; //!< points to the density model
    PhysicsBase *pphy; //!< points to the model of physics
    float *btot,*bpol,*netot; //!< points to the total brightness, polarized brightness and total density images

  private:
    bool *pisrunning;
    bool neonly,quiet,frontinteg;
    unsigned int nbpix,chunksize,lastchunkremain,nbchunk;

  public:
    Camera camera;
    Sun csun;
    Clos los;
    ModelPosition modelposition; //!< model positioning

    friend class PhysicsBase;
    friend class PhysicsThomson;
    friend class PhysicsUV;
    friend class PhysicsFCor;
  
    Scene();
    ~Scene()
    {
        delete pmod,pphy;
    }
    //! Compute the image with multi-threading for each pixel
    void computeImagebyRay(float *btot,float *bpol,float *netot,const unsigned int nbthread);
    //! Compute the image with multi-threading by chunk of the image
    void computeImagebyChunk(float *btot,float *bpol,float *netot,const unsigned int nbthread,const unsigned int nbchunk);
    void setabs(Cbasis abs)
    {
        this->abs=abs;
    };
    void setobs(Cbasis obs)
    {
        this->obs=obs;
    };

    //! Set the density model to be used
    //! \param modelid the id of the model
    //! \param *pmodparam points to the vector of parameters
    void setDensityModel(int modelid,float *pmodparam)
    {
        delete pmod;
        pmod=modelselect(modelid);
        pmod->initParam(pmodparam);
        pphy->setParentScene(this);
    };

    void setNeonly(int neonly)
    {
        this->neonly=(bool)neonly;
    };
    void setQuiet(int quiet)
    {
        this->quiet=(bool)quiet;
    };
    void setFrontInteg(int frontinteg)
    {
        this->frontinteg=(bool)frontinteg;
    };
    void setPhysics(PhysicsType phytype);
    string getPhysics()
    {
        return pphy->getPhysics();
    }
    void setPhysicsParam(float *phyparam) {pphy->setParam(phyparam);}
    void printPhysicsParam() {pphy->printParam();}

private:


//! Thomson scattering integration along a LOS
    inline void losintegThomson(const unsigned int &i,const unsigned int &j)
    {

        Cvec vlosobs=camera.ij2los(float(i),float(j));
        Cvec vlosabs=obs.ui * vlosobs;

        // -- compute rho : dist LOS - Sun cntr: impact parameter
        float rho=psldist(obs.o,vlosabs,abs.o);

        // -- origin of the LOS: observer by default
        Cvec qlos;
        if (frontinteg)
        {
            qlos=obs.o;
        }
        else
        {
            qlos=orthoproj(obs.o,vlosabs,abs.o);
        }

        // -- init LOS integration parameters
        Cvec vlosstep=vlosabs * los.ds;
        Cvec vs=qlos + vlosabs * los.sstart;

        unsigned int pos=i+j*camera.ccd.sxpix;

        for (unsigned int k=0;k<los.nbp;k++,vs+=vlosstep)
        {
            // -- dist Sun cntr - LOS current pos
            float r=vs.norm();

            // -- no integration within and behind the Sun
            if ( (r <= 1.001) || (rho <= 1.001 && (vs-obs.o).mag() > sqrt(obs.o.mag()*obs.o.mag()-rho*rho))) continue;

            // -- Call density here
            float ner=pmod->Density(ChangetoDensityCoord(modelposition,vs));

            if (ner <= 1e-1) continue;

            netot[pos]+=ner;
            if (neonly) continue;

            float btotcoeff,bpolcoeff;
            csun.getThomsonCoeff(r,rho,btotcoeff,bpolcoeff);

            btot[pos] +=ner*btotcoeff;
            bpol[pos] +=ner*bpolcoeff;

        }

        // multiply by the integral constant factor
        btot[pos] *=csun.constfactor*los.ds;
        bpol[pos] *=csun.constfactor*los.ds;
        netot[pos] *=RSUN_CM*los.ds;

    };



//! Generic LOS integration: the physics has to be initialized beforehand
    inline void losinteg(const unsigned int &i,const unsigned int &j)
    {

        Cvec vlosobs=camera.ij2los(float(i),float(j));
        Cvec vlosabs=obs.ui * vlosobs;

        float btout,bpout,neout;
        bool flagnull=0;

        // -- compute rho : dist LOS - Sun cntr: impact parameter
        float rho=psldist(obs.o,vlosabs,abs.o);

        // -- origin of the LOS: observer by default
        Cvec qlos;
        if (frontinteg)
        {
            qlos=obs.o;
        }
        else
        {
            qlos=orthoproj(obs.o,vlosabs,abs.o);
        }

        // -- init LOS integration parameters
        Cvec vlosstep=vlosabs * los.ds;
        Cvec vs=qlos + vlosabs * los.sstart;

        unsigned int pos=i+j*camera.ccd.sxpix;

        for (unsigned int k=0;k<los.nbp;k++,vs+=vlosstep)
        {
            // -- dist Sun cntr - LOS current pos
            float r=vs.norm();

            // -- no integration within and behind the Sun
           // if ( (r <= 1.001) || (rho <= 1.001 && (vs-obs.o).mag() > sqrt(obs.o.mag()*obs.o.mag()-rho*rho))) continue;

            // -- Compute radiation of the volume element depending on chosen physics
            //    Also call the density in here.
            flagnull=pphy->computeRadiation(vs,r,rho,btout,bpout,neout);

            if (flagnull) continue;

            netot[pos]+=neout;
            if (neonly) continue;

            btot[pos] +=btout;
            bpol[pos] +=bpout;
        }

        // ---- multiply by the integral constant factor, depending on the physics
        float btf,bpf,nef;
        pphy->getConstFactors(btf,bpf,nef);
        btot[pos]  *=btf;
        bpol[pos]  *=bpf;
        netot[pos] *=nef;
    }

    //! Multi-threaded Thomson scattering integration along a LOS, for the "by Ray" integration
    void losintegray(const unsigned int &i,const unsigned int &j,const unsigned int &threadid);
    //! Multi-threaded Thomson scattering integration along a LOS, for the "by Chunk" integration
    void losintegchunk(const unsigned int &i,const unsigned int &threadid);

};





#endif