/*
 * Bunch.cc
 *
 *  Created on: May 23, 2011
 *      Author: jakob
 */

#include <math.h>
#include "Bunch.h"
#include "random.h"
#include "Vector3D.h"
#include "Element.h"

namespace vhc {

Bunch::Bunch(const Particle& referenceParticle, int quantity, int lambda, double standardDeviation, double length, double targetEmittance, double A12, double A22):
		Beam(referenceParticle, quantity, lambda),
		standardDeviation(standardDeviation),
		length(length),
		targetEmittance(targetEmittance),
		A12(A12),
		A22(A22) {}

Bunch::~Bunch() {

}

void Bunch::initializeParticles() {
	create(getLength() * referenceParticle.getVelocity().norm());
}

void Bunch::create(double dt) {
	random::seed(0);
	double debit = quantity / length;
	double fraction(debit * dt); // debit "vrai", mais a priori non entier
	int number(fraction);      // partie entière
	fraction -= number;        // partie fractionnaire
	// on ajoute 1 au  hasard, proportionnellement à la partie fractionnaire :
	if ( random::uniform(0.0, 1.0) < fraction ) ++number;

	//calculation des proprietes des particules
	for (int i = 0; i < quantity / lambda; ++i) {
		Particle* particle = getReferenceParticle().clone();

		double sa[2]; // = [r, z]
		double va[2]; // = [vr, vz]
		//calculation de r, vr, z et vz
		for (int i = 0; i < 2; ++i) {
			double A11 = getA11();
			double A12 = getA12();
			double A22 = getA22();

			double theta = atan(2 * A12 / (A11 - A22));
			double a = A11 * cos(theta) * cos(theta) + 2 * A12 * cos(theta) * sin(theta) + A22 * sin(theta) * sin(theta);
			double b = A11 * sin(theta) * sin(theta) + 2 * A12 * cos(theta) * sin(theta) + A22 * cos(theta) * cos(theta);

			double x = random::gaussian(0, sqrt(getTargetEmittance() / a));
			double y = random::gaussian(0, sqrt(getTargetEmittance() / b));

			sa[i] = cos(theta) * x + sin(theta) * y;
			va[i] = -sin(theta) * x + cos(theta) * y;
		}

		double r = sa[0];
		double vr = va[0];

		double z = sa[1];
		double vz = va[1];

		double norm = random::gaussian(getReferenceParticle().getVelocity().norm(), getStandardDeviation());

		Vector3D u = getReferenceParticle().getElement()->getHorizontalAt(getReferenceParticle().getPosition());
		Vector3D v = vr * u  + vz * Vector3D::k + sqrt(norm * norm - vr * vr - vz * vz) * u.cross(Vector3D::k);

		//composante s de position
		double s = random::gaussian(0, getLength());

		//ecart de position
		Vector3D dp = r * u + z * Vector3D::k + s * u.cross(Vector3D::k);

		particle->setMass(getReferenceParticle().getMass() * getLambda());
		particle->setCharge(getReferenceParticle().getCharge() * getLambda());
		particle->translate(dp);

		add(particle);
	}

}

Bunch* Bunch::clone() const {
	return new Bunch(*this);
}


double Bunch::getA11() const {
	return (1 + getA12() * getA12() ) / getA22();
}

double Bunch::getA12() const {
	return A12;
}

double Bunch::getA22() const {
	return A22;
}

double Bunch::getStandardDeviation() const {
	return standardDeviation;
}

double Bunch::getLength() const {
	return length;
}

double Bunch::getTargetEmittance() const {
	return targetEmittance;
}

}