gistfile1.cpp

/*
 *  photon.cpp
 *  photonSplit
 *
 *  Created by Diane Nguyen on 4/25/13.
 *  Copyright 2013 __MyCompanyName__. All rights reserved.
 *
 */
#include "photon.h"
#include <cstdlib>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <iostream>
#include <vector>
using namespace std;

Photon::Photon(){
  random=0;
  cartesienZ=0;
	cosineW=0;
	weight=1.0;
	
}

Photon::Photon(double r, double z, double w, double wt){
	random=r;
	cartesienZ=z;
	cosineW=w;
	weight=wt;
}

double Photon:: getRandom() const{
	return random;
}

double Photon::getCartesienZ() const{
	return cartesienZ;
}

double Photon:: getCosineW() const{
	return cosineW;
}

double Photon:: getWeight() const{
	return weight;
}

void Photon::setRandom(double r){
	random=r;
}

void Photon::setCartesienZ(double z){
	cartesienZ=z;
}

void Photon::setCosineW(double w){
	cosineW=w;
}

void Photon::setWeight(double w){
	weight=w;
}

double Photon::lancement(double xSectionScat, double xSectionPE){
	double z;
	cartesienZ=-log(1-random)/(xSectionScat+xSectionPE);
	return z;
	
}

vector<Photon> Photon:: splitting(int n){
	vector<Photon>split(n);
	double r;
	//srand(time(NULL));
	for(int i=0;i<split.size();i++){
		r=(double)rand()/RAND_MAX;
		//cout<< "random  "<< r<< endl;
		split[i].setRandom(r);
		split[i].setCartesienZ(cartesienZ);
		split[i].setCosineW(cosineW);
		split[i].setWeight(weight/n);
	}
	return split;
}
	   
double Photon::routine(int n, double xSectionScat, double xSectionPE, double thickness){
	bool termination=false;
	double traversedDistance;
	double proba=0.0;
	double z;
	double w;
	double r;
	while(this->getCartesienZ()>0 && !termination){
		if(this->cartesienZ>thickness/2){
			vector<Photon>split= this->splitting(n);
			for(int i=0; i<split.size(); i++){
				bool reaction=true;
				while(split[i].getCartesienZ()<=thickness && reaction){
					//Roulette russe
					if(split[i].getCartesienZ()<=thickness/2){
						split[i].setWeight(split[i].getWeight()/exp(-(xSectionPE+xSectionScat)*split[i].getCartesienZ()));
						split[i].routine(n,xSectionScat,xSectionPE,thickness);
					}
					//PE ou collision elastique?
					r=(double) rand()/RAND_MAX;
					if(r<xSectionScat/(xSectionPE+xSectionScat)){
						traversedDistance=-log(1-split[i].getRandom())/(xSectionScat+xSectionPE);
						w=2*split[i].getRandom()-1;
						z=cartesienZ+w*traversedDistance;
						split[i].setCartesienZ(z);
						split[i].setCosineW(w);
						split[i].setRandom((double) rand()/RAND_MAX);
					}
					else{
						reaction=false;
						termination=true;
					}

				}
				if(split[i].getCartesienZ()>thickness){
					termination=true;
					proba=proba+split[i].getWeight();
				}
			}
		}
		else{
			if(this->getCartesienZ()<thickness/2){
				r=(double) rand()/RAND_MAX;
				if(r<xSectionScat/(xSectionPE+xSectionScat)){
					traversedDistance=-log(1-this->random)/(xSectionScat+xSectionPE);
					w=2*this->random-1;
					z=cartesienZ+w*traversedDistance;
					this->setCartesienZ(z);
					this->setCosineW(w);
					this->setRandom((double) rand()/RAND_MAX);
				}
				else{
					termination=true;
				}
			}
		}
		}
		return proba;
}