/***************************************************************************
 *            losinteg.h
 *
 * $Id: losinteg.h,v 1.1 2009/02/09 20:57:18 thernis Exp $
 ****************************************************************************/


// -- compute pos of ortho proj Sun cntr on the LOS
calcqlos(pofinteg,frontinteg,obs,vlosabs,abs,poi,qlos);

// if (pofinteg == 0 && frontinteg == 0) qlos=orthoproj(obs.o,vlosabs,abs.o); else {
// 	if (pofinteg != 0) {

	// get normal vector to POI: defined as the z axis of POI coord sys
// 	Cvec poinormal=poi.ui.column(3);
	// compute intersection LOS - POI (see 25Aug2005 notes)
// 	qlos=obs.o - vlosabs * (pscal(obs.o,poinormal)/pscal(vlosabs,poinormal));
	/* Integration in the Earth POS (see 23Aug2005 notes)
	float tmpnorm=obs.o.norm();
	qlos=obs.o - vlosabs * ((tmpnorm*tmpnorm) / pscal(obs.o,vlosabs));
	*/

// 	}
// 	if (frontinteg !=0 ) {
// 		// ---- qlos is the observer
// 		qlos=obs.o;
// 	}
// }

// -- initialize
*posbtot=0;
*posbpol=0;
*posne=0;

// -- current position on the LOS
s=los.sstart;
vlosstep=vlosabs * los.ds;
vs=qlos + vlosabs * s;
for(k=0;k<los.nbp;k++,s+=los.ds,vs+=vlosstep) {

  // -- dist Sun cntr - LOS current pos
  r=vs.norm();
  
  // -- integration only in front of the Sun
  if (rho > LIMBLIMIT || (r > LIMBLIMIT && vs.v[2] > 0.)) {
    /******************************************/
    /*            MODEL CALL HERE             */
    /******************************************/
    // -- compute Ne    
    ner=pmod->Density(nps.u * (vs + nps.o), pmodparam,temperature);
    if (ner <= 1e-1) continue;
    
 
    // total density
    *posne +=ner;
 
	// -- skip brightness calculation if requested
	if (flagneonly > 0) continue; 
		
    /******************************************/
    // -- compute geometrical functions
    sinomega=RSUN/r;
    sinsquareomega=sinomega*sinomega;
    cossquareomega=1-sinsquareomega;
    cosomega=sqrt(cossquareomega);
    
    logterm=log((1.+sinomega)/cosomega);
    
    a=cosomega*sinsquareomega;
    b=-1./8.*(1.-3.*sinsquareomega-cossquareomega*((1.+3.*sinsquareomega)/sinomega)*logterm);
    
    c=(4./3.)-cosomega-(cosomega*cossquareomega)/3.;
    d=(1./8.)*(5.+sinsquareomega-cossquareomega*((5.-sinsquareomega)/sinomega)*logterm);
    
    // ---- sum in the pixel
    rhooverr=rho/r;
    // the polarized brightness
    polterm=(a+u*(b-a))*rhooverr*rhooverr;
    *posbpol +=ner*polterm;
    // the total brightness
    *posbtot +=ner*(2*(c+u*(d-c))-polterm);
    
    //if (i == iii && j == jjj) {
    //  *posbtot +=1e9;
    //}
    //	  if (i == 0 && j == 0 && (k == 0 || k == 50)) {
    //if (i == 10 && j == 0) {
    //cout << "sinomega " << sinomega << endl;
    //cout << "cosomega " << cosomega << endl;
    //cout << "a " << a << endl;
    //cout << "b " << b << endl;
    //cout << "c " << c << endl;
    //cout << "d " << d << endl;
    //cout << "sin^2(chi) " << rhooverr*rhooverr << endl;
    //cout << "rho " << rho << endl;
    //cout << "r " << r << endl;
    //cout << "log term " << logterm << endl;
    //cout << "c + u*(d-c)  " << tanterm << endl;
    //cout << "Bt " << tanterm*constfactor*los.ds << endl;
    //cout << "pB " << polterm*constfactor*los.ds << endl;
    //cout << "Br " << (tanterm-polterm)*constfactor*los.ds << endl;
    //cout << "B " << (2*tanterm-polterm)*constfactor*los.ds << endl;
    //cout << "constfactor " << constfactor << endl;
    //cout << "constfactor*ds " << constfactor*los.ds << endl;
    //cout << "u " << u << endl;
    //cout << "ner " << ner << endl;
    //cout << "vs " << vs << endl;
    //cout << "qlos " << qlos << endl;
    //cout << "vlosabs " << vlosabs << endl;
    //cout << "s " << s << endl;
    //cout << "los.ds " << los.ds << endl;
    //cout << "vvv " << vvv << endl;
    //cout << "xx " << xx << endl;
    //cout << "yy " << yy << endl;
    //cout << "xstep " << xstep << endl;
    //cout << "ystep " << ystep << endl;
    //cout << "----------" << endl;
    //}
    
    
  }
}
// multiply by the integral constant factor
*posbtot *=constfactor*los.ds;
*posbpol *=constfactor*los.ds;
*posne *=RSUN_CM*los.ds;



/*
* $Log: losinteg.h,v $
* Revision 1.1  2009/02/09 20:57:18  thernis
* Put some old code here in case I need that one day...
*
* Revision 1.5  2008/09/30 18:33:28  thernis
* Replace include of losinteg.h and losintegadaptstep.h by inline functions.
*
* Revision 1.4  2008/09/23 14:08:12  thernis
* - add new testing models, - implement integration in front of the instrument
*
* Revision 1.3  2008/06/27 18:27:32  thernis
* Fix sing test to integrate only in front of the Sun limb limits, not behind
*
* Revision 1.2  2007/05/14 17:19:40  thernis
* Add CVS id and log in all files
*
*/