#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
#include "/home/gcm/mbenoit/root/include/TMatrixD.h"
#include "/home/gcm/mbenoit/root/include/TMatrixT.h"
#include "/home/gcm/mbenoit/root/include/TH2.h"
#include "/home/gcm/mbenoit/root/include/TApplication.h"
#include "/home/gcm/mbenoit/root/include/TFile.h"









using namespace std;







void initialise(double pix, double gap,double Voltage,int precision,TH2D* V){
	//double per=pix+gap;
	double pas=3*(pix+gap)/precision;
	TMatrixD X(precision,precision);
	TMatrixD Y(precision,precision);
	for(int i=0;i<precision;i++){
		for (int j=0;j<precision;j++){
			X(i,j)=i*pas-3.0*(pix+gap)/2;
			Y(i,j)=j*pas-3.0*(pix+gap)/2;
			};
			};
	
	for(int i=0;i<precision;i++){
		for (int j=0;j<precision;j++){
			//cout << '*' << endl;
			V->SetBinContent(i,j,Voltage);
			/* if (fmod(fabs(X(i,j)+1.5*per),per)<(gap/2+pix) && fmod(fabs(X(i,j)+1.5*per),per)>(gap/2) &&
			fmod(fabs(Y(i,j)+1.5*per),per)<(gap/2+pix) && fmod(fabs(Y(i,j)+1.5*per),per)>gap/2){
				V->SetBinContent(i,j,Voltage);
				}
			else if(Y(i,j)>=X(i,j) && Y(i,j)<=-X(i,j)){
				V->SetBinContent(i,j,(((1+pix/gap) + 2.0*X(i,j)/gap)*Voltage));	
				}
			else if(Y(i,j)>=X(i,j) && Y(i,j)>=-X(i,j)){
				V->SetBinContent(i,j,(((1+pix/gap) - 2.0*Y(i,j)/gap)*Voltage));
				}
			else if(Y(i,j)<=X(i,j) && Y(i,j)<=-X(i,j) ){
				V->SetBinContent(i,j,(((1+pix/gap) + 2.0*Y(i,j)/gap)*Voltage));
				}
			else if(Y(i,j)<=X(i,j) && Y(i,j)>=-X(i,j)){
				V->SetBinContent(i,j,(((1+pix/gap) - 2.0*X(i,j)/gap)*Voltage));
				}; */
			};};
			
			};
			
			
void MasquePix(double pix, double gap,int precision,TMatrixD X, TMatrixD Y,TH2D* V){
		double per=pix+gap;
		for(int i=0;i<precision;i++){
			for (int j=0;j<precision;j++){
				if (fmod(fabs(X(i,j)),per)<pix/2 && fmod(fabs(Y(i,j)),per)<pix/2){
					V->SetBinContent(i,j,0);
					}
				else{
					};
				};};
		};
		
void MasqueGap(double pix, double gap,int precision,TMatrixD X, TMatrixD Y,TH2D* V){
		double per=pix+gap;
		for(int i=0;i<precision;i++){
			for (int j=0;j<precision;j++){
				if (fmod(fabs(X(i,j)),per)<pix/2 && fmod(fabs(Y(i,j)),per)<pix/2){
				
					}
				else{
					V->SetBinContent(i,j,0);
					};
				};};
				
		};
void ConditionV(double pix, double gap,int precision,double Voltage,TMatrixD X, TMatrixD Y,TH2D* V){
		double per=pix+gap;
		for(int i=0;i<precision;i++){
			for (int j=0;j<precision;j++){
				if(V->GetBinContent(i,j)>Voltage){
					V->SetBinContent(i,j,Voltage);
					};
				if(fmod(fabs(X(i,j)+1.5*per),per)<(gap/2+pix) && fmod(fabs(X(i,j)+1.5*per),per)>(gap/2) &&
			fmod(fabs(Y(i,j)+1.5*per),per)<(gap/2+pix) && fmod(fabs(Y(i,j)+1.5*per),per)>gap/2){
					V->SetBinContent(i,j,Voltage);
					};
			};};
				
		};
		
				

void boucle(double pix,double gap,double Voltage, int precision, TH2D* Vprec,TMatrixD X, TMatrixD Y){
	double L=1.0,pi=3.1416,per=pix+gap;
	TH2D *E=new TH2D(*Vprec);
	E->FFT(E,"RE DCT_0 M");
	E->Scale(1.0/precision);
	TH2D *Epix=new TH2D(*E);
	TH2D *Egap=new TH2D(*E);
	for(int i=0;i<precision;i++){
		for (int j=0;j<precision;j++){
			if(i==0 && j==0){Epix->SetBinContent(i,j,E->GetBinContent(i,j)*(1.0/L));
					 Egap->SetBinContent(i,j,E->GetBinContent(i,j)*(1.0/L));}
			else{
		
		
		
Epix->SetBinContent(i,j,E->GetBinContent(i,j)*((1+exp(-2*(2*pi*sqrt(double(i*i+j*j))/(precision-1))*L))/(1-exp(-2*(2*pi*sqrt(double(i*i+j*j))/(precision-1))*L)))*(2*pi*sqrt(double(i*i+j*j))/(precision-1)));
			Egap->SetBinContent(i,j,E->GetBinContent(i,j)*((2*pi*sqrt(double(i*i+j*j)))/(precision-1)));

			};
			};
			};
	
	Epix->FFT(Epix,"RE DCT_0 M");
	Epix->Scale(1.0/precision);
	Egap->FFT(Egap,"RE DCT_0 M");
	Egap->Scale(1.0/precision);
	TH2D *Etotal = new TH2D("Etotal","Etotal",precision,-1.5*per,1.5*per,precision,-1.5*per,1.5*per);
	MasquePix(pix,gap,precision,X,Y,Egap);
	MasqueGap(pix,gap,precision,X,Y,Epix);
	Etotal->Add(Epix,Egap,1,1.0/14);
	Etotal->FFT(Etotal,"RE DCT_0 M");
	Etotal->Scale(1.0/precision);
	TH2D *Vf=new TH2D(*Etotal);
	for(int i=0;i<precision;i++){
		for (int j=0;j<precision;j++){
			if(i==0 && j==0){Vf->SetBinContent(i,j,Etotal->GetBinContent(i,j)*L);}
			else{				
	Vf->SetBinContent(i,j,Etotal->GetBinContent(i,j)/(((1+exp(-2*(2*pi*sqrt(double(i*i+j*j))/(precision-1))*L))/(1-exp(-2*(2*pi*sqrt(double(i*i+j*j))/(precision-1))*L)))*(pi*sqrt(double(i*i+j*j))*L/(precision-1))));
			};
			};
			};
	
				
	Vf->FFT(Vf,"RE DCT_0 M");
	Vf->Scale(1.0/precision);
	ConditionV(pix,gap,precision,Voltage,X,Y,Vf);
	Vf->Scale(0.1);
	Vf->Add(Vprec,0.9);
	for(int i=0;i<precision;i++){
		for (int j=0;j<precision;j++){
			Vprec->SetBinContent(i,j,Vf->GetBinContent(i,j));
			};};
	delete Vf;
	delete Etotal;
	delete E;
	delete Epix;
	delete Egap;

	};
	

int solvercompile(){
	TFile f("histofinal.root","RECREATE");
	
	int precision=1000,b=0;
	double gap=100e-6,pix=125e-6,param=10000,per=pix+gap,pas=3*per/precision;
	double Voltage=1100;
	TMatrixD X(precision,precision);
	TMatrixD Y(precision,precision);
	for(int i=0;i<precision;i++){
		for (int j=0;j<precision;j++){
			X(i,j)=i*pas-3*(pix+gap)/2;
			Y(i,j)=j*pas-3*(pix+gap)/2;
			};
			};
	TH2D V("V","V",precision,-1.5*per,1.5*per,precision,-1.5*per,1.5*per);
	
	initialise(pix,gap,Voltage,precision,&V);
 	while(param>1e-5){
	TH2D Vp(V);
	boucle(pix,gap,Voltage,precision,&V,X,Y);
	param=fabs(V.Integral() - Vp.Integral());
	cout << param << endl;
	b++;
	};
	V.SetMinimum(0); 
	V.SetMaximum(1100);
	V.Draw("LEGO");
	V.Write();
	

	return 0;	
	
	}