#include <iostream>
#include <fstream>
using namespace std;

#include <math.h>
#include "beam.h"
#include "starlightconstants.h"
#include "inputparameters.h"
#include "bessel.h"
Beam::Beam(int Zin, int Ain, double bdeuteron, int in_or_co, Inputparameters& input) 
: Nucleus(Zin,Ain,bdeuteron,in_or_co)
{

	//Setting needed inputparameters to protected variables
	BeamInputGamma_em=input.getgamma_em();

//just a lorentz boost for the time being, maybe make a separate function for 
//this and the constructor pretty harmless?  
//Otherwise, we will need to insert momenta
//so for the moment it does nothing, 
//but has the boost in place.do it as an array?
//Which gamma? Egamma or gamma_em, I say Egamma? 
//Either energy of the photon or the lorentz value from the beam.

//If done right, this might allow for linear accelerators?

/*
	double betax,betay,betaz;
	double Ecm,W;
	double px0,py0,pz0;
	Egamma = gamma;
	Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
	betax = -px0/Ecm;
	betay = -py0/Ecm;
	betaz = -pz0/Ecm;
*/

}

double Beam::nofe(double impactparameter)
{
//Function for the calculation of the "photon density".
//nofe=numberofgammas/(energy*area)
//Assume beta=1.0 and gamma>>1, i.e. neglect the (1/gamma^2)*K0(x) term
    double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.;
    
    X = (impactparameter*Egamma)/(BeamInputGamma_em*StarlightConstants::hbarc);
    if( X <= 0.0) cout<<"In nofe, X= "<<X<<endl;
    
    factor1 = (double(getZin()*getZin())*StarlightConstants::alpha)/
    (StarlightConstants::pi*StarlightConstants::pi);
    
    factor2 = 1./(Egamma*impactparameter*impactparameter);
    factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X));
    nofex    = factor1*factor2*factor3;


    return nofex;



}

Beam::~Beam()
{
}
