From 9c61b7a85abfa2ed5668f7932529eade0f37113e Mon Sep 17 00:00:00 2001
From: Peter DeRosa <54421997+MagnusSnowleopard@users.noreply.github.com>
Date: Mon, 6 Jun 2022 12:35:44 -0400
Subject: [PATCH] Add files via upload

---
 1062.csv |   4 +
 AD.cxx   | 750 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Cleb.cxx | 740 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Cleb.h   | 653 ++++++++++++++++++++++++++++++++++++++++++++++++
 Coeff.h  |  71 ++++++
 GUI.h    | 422 +++++++++++++++++++++++++++++++
 QDK.h    | 206 +++++++++++++++
 Racah.h  | 303 ++++++++++++++++++++++
 ad.txt   |   7 +
 9 files changed, 3156 insertions(+)
 create mode 100644 1062.csv
 create mode 100644 AD.cxx
 create mode 100644 Cleb.cxx
 create mode 100644 Cleb.h
 create mode 100644 Coeff.h
 create mode 100644 GUI.h
 create mode 100644 QDK.h
 create mode 100644 Racah.h
 create mode 100644 ad.txt

diff --git a/1062.csv b/1062.csv
new file mode 100644
index 0000000..98dd2da
--- /dev/null
+++ b/1062.csv
@@ -0,0 +1,4 @@
+angle,Y,Yerr
+45,129,10
+90,110,10
+135,129,10
diff --git a/AD.cxx b/AD.cxx
new file mode 100644
index 0000000..368144d
--- /dev/null
+++ b/AD.cxx
@@ -0,0 +1,750 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+
+#include "GUI.h"
+#include "Coeff.h"
+#include "QDK.h"
+#include "Racah.h"
+#include "Cleb.h"
+
+
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+//For more information about angular distributions, read Rose and Brinks, and Frank Moore (1988). 
+
+using namespace std;
+void menu(){
+
+	std::cout <<"==========================================\n";
+	std::cout<< "\t\t Menu Options \t \n";
+	std::cout<< "0 - Reading and Instructions \n"; 
+	std::cout<< "1 - Plot Chi Sqr \n";
+	std::cout<< "2 - Plot Ang.Dis Fit \n"; 
+	std::cout<< "3 - Clear J1 J2 Memory \n";
+	std::cout<< "4 - Exit \n";
+	std::cout<< "5 - Legfit to Data \n"; 
+	std::cout<< "6 - Generate Data \n"; 
+	std::cout<<"=========================================== \n";} 
+void Readme(){
+	std::cout<< " The program calculates Chi-Squared values \n";
+	std::cout<< " from experimental angular distributions as \n";
+	std::cout<< " a function of multipole ratios using the \n"; 
+	std::cout<< " theoretical angular distribution formulae \n"; 
+	std::cout<< " in rose and brink \n";
+	std::cout<< " \n"; 
+	std::cout<< " Follow the prompt in order to correctly display \n";
+	std::cout<< " \n";
+	std::cout<< " To close the gui, press any button. \n";
+	std::cout<< " To zoom in, left click then drag and let go. To unzoom\n";
+	std::cout<< " press the space bar. To draw, right click\n";
+	std::cout<< "                                          \n";
+	std::cout<< " ad.txt is generated with geometric stats.\n";
+	std::cout<< " Then at the end it generates CG, W3J & W6J\n";
+	
+}
+#ifndef __CINT__ 
+
+int main(int argc,char **argv){
+	HistoGUI gui;
+	//	menu();
+	double j1, j2;
+
+
+	//detector data input
+	double detradius, targetdistance, detthickness = -1;
+	double Energy;
+	//input .csv file
+	double Sigma; 	
+	double Feeding; 
+
+//--------------------------------
+
+	detradius = 3.;
+	targetdistance = 4.;
+	detthickness = 5.;
+	
+	Energy = 1062.;
+	Sigma = .1; 
+	Feeding = .1;
+
+	j1 = 2.; 
+	j2 = 1.;
+//--------------------------------
+	
+//TEST MODE? y : 1 || n = 0
+
+	int test = 1; 
+
+	int det_param_token = -1;
+	int gamma_energy_token = -1;
+	int ang_file_token = -1; 
+	int sigma_token = -1; 
+	int feeding_token = -1;
+	int j1j2token = -1; 
+ if(test == 1){
+
+	 det_param_token = 1;
+	 gamma_energy_token = 1;
+	 ang_file_token = -1; 
+	 sigma_token = 1; 
+	 feeding_token = 1;
+	 j1j2token = 1; 
+
+}
+/*
+	printf("Input detector radius, target distance, and detector thickness : ");
+	scanf("%lf,%lf,%lf",&detradius,&targetdistance,&detthickness);
+
+
+	if(detradius > 0 && targetdistance > 0 && detthickness > 0){
+		det_param_token = 1;
+	}else{
+
+		do{
+			if(detradius< 0|| targetdistance < 0|| detthickness < 0){
+				printf("Negative numbers are not allowed!\nRe-enter radius, distance, and thickness : ");
+				scanf("%lf,%lf,%lf", &detradius,&targetdistance,&detthickness);
+			}
+		}while(detradius< 0|| targetdistance < 0|| detthickness < 0);}
+
+	if(detradius > 0 && targetdistance > 0 && detthickness > 0){
+		det_param_token = 1; printf(" --- Detector characteristics loaded --- \n");
+	}
+
+	printf("Enter the gamma-ray energy : ");
+	scanf("%lf", &Energy);
+	if(Energy > 0){
+		gamma_energy_token = 1;
+	}else{
+		do{
+			if(Energy <= 0){		
+				printf("Negative numbers are not allowed!\nRe-enter gamma-ray energy : ");
+			scanf("%lf", &Energy);}
+
+		}while(Energy <= 0);}
+
+	if(Energy > 0){
+		gamma_energy_token = 1; printf(" --- Gamma Energy loaded --- \n");
+	}*/
+
+	QKN(Energy, detradius, targetdistance, detthickness);
+	//then calucate QD2 and QD4, and replace the 0 with them. 
+	double QD2 = .7;
+	double QD4 = .3;
+
+	//input data file of Angular Data (theta, Yexp, Yerr)
+	string fname;
+	cout<<"Enter the file name : ";
+	cin>>fname;
+
+	vector<vector<string> > content;
+	vector<string> row;
+	string line, word; 
+	
+	vector<string> adata;
+	vector<string> ydata;
+	vector<string> eydata; 
+	int num = 0;
+	fstream file(fname.c_str());
+	if(file.is_open()){
+		while(getline(file,line)){
+			row.clear();
+	
+			stringstream str(line);
+
+			while(getline(str, word, ','))
+			
+				row.push_back(word);
+				content.push_back(row);
+				adata.push_back(row[0]);
+				ydata.push_back(row[1]);
+				eydata.push_back(row[2]);			
+				 
+		}
+		cout<< " --- Angular Data File Loaded --- \n";
+	}else cout<< "Could not open the file\n";
+	
+	for(int i = 0; i< adata.size();i++){
+//		cout << adata[i] << "\n";
+//		cout << ydata[i] << "\n";
+//		cout << eydata[i] << "\n";
+	}
+	
+	// Legendre Polinomial fit using adata, ydata, eydata;
+
+	vector<int> int_adata(adata.size());
+std::transform(adata.begin()+1, adata.end(), std::back_inserter(int_adata), StrToInt);
+
+	vector<double> double_angle(int_adata.size() -4);
+	for(int i = 4; i<  int_adata.size(); i++){
+	//	printf("angle = %d\n",int_adata[i]);
+		double_angle.push_back((double)int_adata[i]*3.14159/180);	
+	}
+
+	vector<int> int_ydata(ydata.size());
+std::transform(ydata.begin()+1, ydata.end(), std::back_inserter(int_ydata), StrToInt);
+
+	for(int i = 4; i<  int_ydata.size(); i++){
+	//	printf("Ydata = %d\n",int_ydata[i]);
+	}
+
+	vector<double> double_ydata(int_ydata.size() -4);
+	for(int i = 4; i<  int_ydata.size(); i++){
+	//	printf("angle = %d\n",int_adata[i]);
+		double_ydata.push_back((double)int_ydata[i]);	
+
+	}
+	
+	vector<int> int_eydata(eydata.size());
+std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), StrToInt);
+
+	for(int i = 4; i<  int_eydata.size(); i++){
+	//	printf("Error data = %d\n",int_eydata[i]);
+	}
+
+	
+	for(int i = 3; i< double_angle.size();i++){
+	//	printf("rad angle = %lf\n",double_angle[i]);
+	
+		//displays the angles inputed as doubles in radians. 
+	}
+	
+	for(int i = 3; i< double_ydata.size();i++){
+	//	printf("y_data = %lf\n",double_ydata[i]);
+	
+		//displays the angles inputed as doubles in radians. 
+	}
+	/*
+	for(int i=0; i<content.size(); i++){
+		for(int j=0; j<content[i].size(); j++){
+			cout<< content[i][j]<<" ";
+		}
+		cout<<"\n";
+	}
+	*/
+	//legendre fitting 
+	// F[xi] = A0 + A2(P2(xi)) + A4(P4(xi))
+	//xi data is angle data in radians. 
+	//yi data is y-intensity data. 
+	
+	//to change the number of angles, just mess with the indexing of the read in and arrays. 
+
+
+	//12 constants. this will build our matrix for gaussian elinmination. 
+	//[a1 a2   a3] [A0] = [a4]
+	//[a5 a6   a7] [A2] = [a8]
+	//[a9 a10 a11] [A4] = [a12]
+	
+	double a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12 = 0.;
+	
+	for(int i = 3; i< double_angle.size(); i++){
+		
+//eq1
+		a1 += 1; //how many data points, not intensional but convientient.  
+		a2 += (1.5 * pow(double_angle[i],2) - 1);
+		a3 += (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a4 += double_ydata[i]; 
+//eq2 
+		a5 += (1.5 * pow(double_angle[i],2) - 1);
+		a6 += (1.5 * pow(double_angle[i],2) - 1)*(1.5 * pow(double_angle[i],2) - 1);
+		a7 += (1.5 * pow(double_angle[i],2) - 1)*(35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a8 += (1.5 * pow(double_angle[i],2) - 1)*double_ydata[i];
+//eq4
+		a9 += (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a10 += (1.5 * pow(double_angle[i],2) - 1)*(35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a11 +=(35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. )* (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a12 += (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. )*double_ydata[i];
+
+		 
+	}
+/*
+	printf("a1 = %lf\n",a1);	
+	printf("a2 = %lf\n",a2);	
+	printf("a3 = %lf\n",a3);	
+	printf("a4 = %lf\n",a4);	
+	printf("a5 = %lf\n",a5);	
+	printf("a6 = %lf\n",a6);	
+	printf("a7 = %lf\n",a7);	
+	printf("a8 = %lf\n",a8);	
+	printf("a9 = %lf\n",a9);	
+	printf("a10 = %lf\n",a10);	
+	printf("a11 = %lf\n",a11);	
+	printf("a12 = %lf\n",a12);	
+*/
+//here number of equations are 3 because its a 3 coeff. fit. 
+	int n = 3; 
+	double mat[3][4] = {{a1,a2,a3,a4},{a5,a6,a7,a8},{a9,a10,a11,a12}};
+//where A0,A2,A4 are stored
+	double residual[n];
+
+	for(int i = 0; i<n;i++){
+		for( int j = i+1; j<n; j++){
+			if(abs(mat[i][i]) < abs(mat[i][j])){
+				for(int k = 0; k<n+1;k++){
+					//swapping matrix
+		//			mat[i][k] = mat[i][k] + mat[j][k];
+		//			mat[j][k] = mat[i][k] - mat[j][k];
+		//			mat[i][k] = mat[i][k] - mat[j][k];
+					swap(mat[i][k],mat[j][k]);
+				}
+
+			}
+		}
+	}
+//performing gaussian elimination.
+
+	for(int i = 0; i< n-1; i++) {
+
+		for(int j = i +1; j<n;j++){
+			
+			float f = mat[j][i]/mat[i][i];
+
+			for(int k = 0; k<n+1; k++){
+				mat[j][k] = mat[j][k] -f*mat[i][k];
+
+			}
+		}
+	}
+	
+//backwards substitution for unknowns.
+
+	for( int i = n-1; i >= 0; i--){
+
+		residual[i] = mat[i][n];
+
+		for( int j = i + 1; j<n;j++){
+
+			if(i != j){
+				residual[i] = residual[i] - mat[i][j] * residual[j];
+			}
+		}
+		residual[i] = residual[i]/mat[i][i]; 
+	}
+	
+
+	cout<<"A0 = "<<residual[0]<<"\n";
+	cout<<"A2 = "<<residual[1]<<"\n";
+	cout<<"A4 = "<<residual[2]<<"\n";
+
+	
+
+
+	/*
+
+	printf("Please enter sigma, sigma = 0 for perfect allignment : ");
+	scanf("%lf", &Sigma);
+	if(Sigma >= 0 ){
+		sigma_token = 1; 
+	}else{
+		do{
+
+			printf("Negative numbers are not allowed!\nRe-enter Sigma : ");
+			scanf("%lf", &Sigma);
+		}while(Sigma <0);}
+
+	if(Sigma >= 0 ){
+		sigma_token = 1; printf(" --- Sigma Loaded --- \n"); 
+	}
+
+	printf("Please enter feeding parameter : ");
+	scanf("%lf", &Feeding);	
+	if(Feeding >= 0 and Feeding < 1){
+		feeding_token = 1;
+	}else{
+		do{
+			printf("Negative numbers are not allowed!\nRe-enter Feeding : ");
+			scanf("%lf", &Feeding);
+		}while(Feeding < 0);}
+
+	if(Feeding >= 0 ){
+		feeding_token = 1; printf(" --- Feeding Loaded --- \n");
+	}
+
+	printf("Please enter J1,J2 : ");
+	scanf("%lf,%lf",&j1,&j2);	
+	if(j1 >= 0 && j2  >= 0 ){
+		j1j2token = 1;
+	}else{
+		do{
+			printf("Negative numbers are not allowed!\nRe-enter J1,J2 : ");
+			scanf("%lf,%lf", &j1,&j2);	
+		}while(j1 < 0 || j2 < 0);}
+
+	if(j1 >= 0 && j2  >= 0 ){
+		j1j2token = 1; printf(" --- J1 & J2 Loaded --- \n");
+	}
+*/
+	
+
+	/*
+	   printf("detector radius = %lf\n",detradius);
+	   printf("detector distance = %lf\n",targetdistance);
+	   printf("detector width = %lf\n",detthickness);
+
+	   printf("Gamma-Ray energy = %lf\n", Energy);
+	   printf("Sigma Distribution = %lf\n", Sigma); 
+	   printf("Feeding param = %lf\n", Feeding); 
+	   printf("J1,J2 = %lf,%lf\n",j1,j2);
+	   */
+/*	
+	double J  = j1 + j2;
+	double m1 = j1;
+	double m2 = j2;
+	double M  = m1 + m2; 
+	double j3 = J; 
+	double m3 = j3;	
+	double j4 = j3;
+	double j5 = j1 + j2 + j3;
+	double j6 = j2 + j3; 
+*/
+// Calculate BK(j1) for perfect allignment
+//
+	double A[6];
+	double js = sqrt(2*j1 +1); 
+	double j12 = 2*j1;
+	int IS = (int)j12 % 2; 
+	
+	double cg12,cg22,cg24,cg14 = 0.;
+	double I1, I2 = 0.;
+	double Bk11, Bk12 = 0.;
+	
+	double cg1, cg2 = 0.;	
+	
+	if(IS == 1){
+		A[0] = j1;
+		A[1] = j1;
+		A[2] = 2.;
+		A[3] = .5;
+		A[4] = -.5;
+		A[5] = 0.;
+	
+		cg12 = CG2(A);
+		
+		A[3] = -.5;
+		A[4] = .5;
+		
+		cg22 = CG2(A);
+
+		A[2] = 4.;
+		
+		cg24 = CG2(A);
+
+		A[3] = .5;
+		A[4] = -.5;
+
+		cg14 = CG2(A);
+
+		Bk11 = (0.5)*js*(pow(-1,I1) * cg12 + pow(-1,I2) * cg22);
+		Bk12 = (0.5)*js*(pow(-1,I1) * cg14 + pow(-1,I2) * cg24);
+	}else if(IS == 0){
+		A[0] = j1; 
+		A[1] = j1; 
+		A[2] = 2.;
+		A[3] = 0.;
+		A[4] = 0.;
+		A[5] = 0.; 
+		
+		cg1 = CG2(A);
+	
+		Bk11 = pow(-1,j1) * js * cg1;
+
+		A[2] = 4.;
+	
+		cg2 = CG2(A);
+		Bk12 = pow(-1,j1) * js * cg2; 	
+		
+	}
+
+//normalize W(m1) 
+//
+
+	       j12 = 2*j1;
+	double j14 = 4*j1;
+	
+	double sigsq = pow(Sigma,2);
+	double sum1 = 0.; 
+
+	double am1,amsq,x,ex = 0.;
+
+	for(int i = 0; i <= j14; i = i + 2){
+		am1 = 0.5*(i - j12);
+		amsq = pow(am1,2);
+		x = - (amsq/(2*sigsq));
+		ex = exp(x);
+
+		sum1 = sum1 + ex;
+	}
+	double cn1 = 1./sum1;	
+
+	double AL0 = 0.;
+	double A0 = j1 - j2; 
+	double A0b = abs(A0);
+	if(A0b > 0){ AL0 = A0b;}
+	
+	if(A0b <= 0){AL0 = 1.;}
+
+	double AL1 = AL0 +1; 
+
+	double am11,amsq1,x1,ex1 = 0.;
+//calculate Bk(j1) for gaussian W(m1) or non zero Sigma
+
+	A[0] = j1; 
+	A[1] = j1; 
+	
+	double sfact = sqrt(2*j1 +1); 
+	double cgg = 0.; 
+	double tTerm = 0.; 
+	double II = 0.;
+	for(int i = 2; i<= 4; i = i + 2){
+		
+		A[2] = i; 
+		Bk11 = 0.;
+		Bk12 = 0.; 
+		for(int m = 0; m <= j14; m = m +2){
+			am11 = 0.5*(m - j12);
+			amsq1 = pow(am11,2);
+			x1 = -(amsq/(2*sigsq));
+			II = j1 - am11; 
+			A[3] = am11; 
+			A[4] = -am11; 
+			cgg = CG2(A);
+			ex1 = exp(x1);
+			tTerm = cn1*ex1*pow(-1,II) * sfact*cgg;
+			if(i == 2) Bk11 = Bk11 * tTerm;
+			if(i == 4) Bk12 = Bk12 * tTerm; 
+		}
+	}
+	
+
+//need to add correction of statistical tensors for cascade feeding. 
+//
+//
+
+
+//Calculate Rk(LL,j1j2)
+//
+	double B[6];
+	
+	double c0,c1,c2,q0,q1,q2,i0,i1,i2,sf0,sf1,sf2, rk01,rk11,rk21, rk02, rk12,rk22 = 0.; 
+	
+	for(int k = 1; k < 4; k = k + 2){
+
+		A[0] = AL0;
+		A[1] = AL0;
+		A[2] = k; 
+		A[3] = 1.;
+		A[4] = -1.;
+		A[5] = 0.;
+		
+		B[0] = j1; 
+		B[1] = j1; 
+		B[2] = AL0; 
+		B[3] = AL0; 
+		B[4] = k; 
+		B[5] = j2; 
+	
+		c0 = CG2(A);
+		q0 = RACAH(B);
+		i0 = 1+j1 -j2-k; 
+
+		sf0 = sqrt((2*j1+2) * (2*AL0 +1) * (2*AL1 +1)); 
+
+		if(k == 2){ rk01 = (pow(-1,i0) * sf0 * c0 * q0);}	
+		if(k == 4){ rk02 = (pow(-1,i0) * sf0 * c0 * q0);}
+
+		A[1] = AL1; 
+		B[3] = AL1; 
+
+		c1 = CG2(A);
+		q1 = RACAH(B);
+		i1 = 1 +j1 - j2 + AL1 - AL0 - k;
+		sf1 = sqrt((2*j1 +1) * (2*AL0 +1) *(2*AL1 +1));
+		 
+		if(k == 2){ rk11 = (pow(-1,i1) * sf1 * c1 * q1);}
+		if(k == 4){ rk12 = (pow(-1,i1) * sf1 * c1 * q1);}
+
+		A[0] = AL1; 
+		B[2] = AL1; 
+		
+		c2 = CG2(A); 
+		q2 = RACAH(B); 
+		i2 = 1 + j1 - j2 - k;
+		sf2 = sqrt((2*j1 +1) * (2*AL1+1) *(2*AL1+1));
+
+		if(k == 2){ rk21 = (pow(-1,i2) * sf2 * c2 * q2);}
+		if(k == 4){ rk22 = (pow(-1,i2) * sf2 * c2 * q2);}
+		
+
+	}
+
+
+
+
+
+	
+	printf("Bk1 = %lf\n",Bk11);
+	printf("Bk2 = %lf\n",Bk12);
+
+	printf("rk01 = %lf\n",rk01); 
+	printf("rk11 = %lf\n",rk11); 
+	printf("rk21 = %lf\n",rk21); 
+	printf("rk02 = %lf\n",rk02); 
+	printf("rk12 = %lf\n",rk12); 
+	printf("rk22 = %lf\n",rk22); 
+	
+
+/*
+
+	double JJ[6] = {j1,j1,2,.5,-.5,0};
+	double racah = RACAH(JJ);	
+	double cg222 = CG2(JJ);
+
+	
+	
+//	double cg = CG(J,M,j1,m1,j2,m2);
+//	double W3J= ThreeJsym(j1,m1,j2,m2,j3,m3);
+//	double W6J = SixJsym(j1,j2,j3,j4,j5,j6);
+
+	double cg_a = CG_a(JJ);
+	double W3J_a= ThreeJsym_a(JJ);
+	double W6J_a = SixJsym_a(JJ);
+	printf(" CG convert : %lf\n",cg222);	
+	printf(" CG         : %lf\n", cg_a);
+	printf(" W3-J 	    : %lf\n", W3J_a);
+	printf(" W6-J       : %lf\n", W6J_a);	
+	printf(" Racah convert  : %lf\n", racah);
+
+*/	
+	// here the Chi-sqs from A2 and A4 needed to be added together. 
+
+	 
+	double delta_min = -3.14159 / 2.;
+	double delta_max = 3.14159 / 2.;
+	double step = 0.001; 
+
+ 	double delta = 0.;
+	double tan_delta =0.;
+	double A0E = residual[0];//NEED FrOM AF FIT
+	double A2E = residual[1];
+	double A4E = residual[2]; 
+	double A2T,a2T = 0.;
+	double A4T,a4T = 0.;
+	double a4E = residual[2]; //NEED FROM AD FIT
+	double rd2T = 0.;
+	double rd4T = 0.;
+	double X22,X24 = 0.;
+	double X2_total = 0.;
+
+	int points = (delta_max - delta_min) / step ; 
+	
+	vector<double> chisqr;
+	vector<double> tdelta;
+
+	for(int i = 0; i< points; i++){
+		tan_delta = i + step + delta_min;
+		delta = tan(tan_delta);
+		rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
+		A2T = QD2*Bk11*rd2T;
+		a2T = A2T/A0E;
+
+		
+		rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2));
+		A4T = QD4*Bk12*rd4T;
+		a4T = A4T/A0E;
+	
+		X22 = pow(abs(A2E- a2T),2)/(abs(a2T));
+		X24 = pow(abs(A4E -a4T),2)/(abs(a4T));
+		X2_total = X22 + X24;
+
+		chisqr.push_back(X2_total);
+		tdelta.push_back(tan_delta);
+
+	}
+
+
+
+
+
+
+	int param;
+
+	param = param_run(det_param_token, gamma_energy_token, ang_file_token, sigma_token, feeding_token, j1j2token);  	
+
+
+	int optnum = -1;
+	menu();
+
+	printf("Input option : ");
+	scanf("%d", &optnum);
+	if(optnum < 0){
+		optnum = -1;
+		printf("Negative numbers are not allowed!\nInput option : ");
+		scanf("%d", &optnum);
+	}
+
+	if(optnum == 0){
+		Readme();
+		optnum = -1;
+
+		printf("Input option : ");
+		scanf("%d", &optnum);
+	} 
+
+	//if all inputs are valid, plot distribution. 
+	if(param == 1 && optnum == 1){
+		optnum = -1;
+
+		//plotting values here
+		std::vector<double> x;
+		std::vector<double> y;
+		
+	        //x will be arctan of mixing ratio. 
+	        //y will be log(X^2); 
+	        
+				
+
+		for(int i = 0; i < 50; i++){
+			x.push_back( i );
+			y.push_back( pow(i, 0.5) );
+		}
+		gui.SetData(x,y);
+		gui.Init();
+		gui.Loop();
+		gui.Close();
+/*
+		printf("Input option : ");
+		scanf("%d", &optnum);
+		if(optnum < 0){		
+			printf("Input option : ");
+			scanf("%d", &optnum);
+		}*/	
+	}
+
+
+
+	return 1;
+
+
+}
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Cleb.cxx b/Cleb.cxx
new file mode 100644
index 0000000..4c3796c
--- /dev/null
+++ b/Cleb.cxx
@@ -0,0 +1,740 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+
+#include "Coeff.h"
+
+using namespace std;
+
+
+
+double CG2(double A[6]){
+
+	double ii[6];
+	double im[6];
+	double ix[9];
+	double nx[5];
+
+	for(int i = 0;i<6;i++){
+		ii[i] = 2.0*A[i];
+	}
+	double Cleb = 0.;
+	double IA = ii[0];
+	double IB = ii[1];
+	double IC = ii[2];
+	double ID = ii[3];
+	double IE = ii[4];
+	double IF = ii[5];
+
+	if(ID + IE - IF != 0){return 0;}
+	double K0 = IA + IB + IC;
+	if((int)K0 % 2 !=0){return 0;}
+
+	double K1 = IA + IB - IC;
+	double K2 = IC + IA - IB;
+	double K3 = IB + IC - IA;
+	double K4 = IA - abs(IB-IC);
+	double K5 = IB - abs(IC-IA);
+	double K6 = IC - abs(IA-IB);
+
+	double karr[6]  = {K1,K2,K3,K4,K5,K6};
+	double kmn=karr[0];
+	double kmx=karr[0];
+	for(int i=0;i<6;i++)
+	{
+		if(kmn>karr[i])
+		{
+			kmn=karr[i];
+		}
+		else if(kmx<karr[i])
+		{
+			kmx = karr[i];
+		}
+	}
+	double K7 = kmn; 
+/*
+	printf("K1 = %lf\n",K1);
+	printf("K2 = %lf\n",K2);
+	printf("K3 = %lf\n",K3);
+	printf("K4 = %lf\n",K4);
+	printf("K5 = %lf\n",K5);
+	printf("K6 = %lf\n",K6);
+	printf("K7 = %lf\n",K7);
+*/
+
+	if(K7 <  0) return 0; 
+
+	for(int i = 0;i<3;i++){
+		if((int)(ii[i]+ii[i+3])%2 !=0)return 0;
+		if(ii[i] < abs(ii[i+3]))return 0; 
+	}
+	double sgnfac =1.0;
+
+	for(int i = 0; i<6;i++){
+		im[i] = ii[i];
+	}
+	double FC2, IT = 0.;
+	//-------------------------------------------------------------------------
+	if(IA >= IB){
+		if(IC >= IB){
+			if(IB != 0){
+				if(IE < 0){
+				//assuming the logic passes;
+					FC2 = IC +1;
+
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+
+
+
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+					printf("Here place 1\n");
+					printf("ix[0] = %lf\n",ix[0]);
+					printf("ix[1] = %lf\n",ix[1]);
+					printf("ix[2] = %lf\n",ix[2]);
+					printf("ix[3] = %lf\n",ix[3]);
+					printf("ix[4] = %lf\n",ix[4]);
+					printf("ix[5] = %lf\n",ix[5]);
+					printf("ix[6] = %lf\n",ix[6]);
+					printf("ix[7] = %lf\n",ix[7]);
+					printf("ix[8] = %lf\n",ix[8]);
+	
+					printf("sqdclg = %lf\n",sqfclg);
+
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+				}else{
+	// if (IE >=0) _>>>
+					ID = -ID;
+					IE = -IE;
+					IF = -IF;
+					if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1)*sgnfac;
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					FC2 = IC + 1; 
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					
+					printf("IT = %lf\n",IT);
+					printf("sqfclg = %lf\n",sqfclg);
+					
+					printf("FC2 = %lf\n",FC2);
+					printf("lgamma(IT) = %lf\n",lgamma(IT+1));
+					printf("log(fC2) = %lf\n",log(FC2));
+
+					printf("IA = %lf\n",IA);
+					printf("IB = %lf\n",IB);
+					printf("IC = %lf\n",IC);
+					printf("ID = %lf\n",ID);
+					printf("IE = %lf\n",IE);
+					printf("IF = %lf\n",IF);
+
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+						printf("sqfclg = %lf\n",sqfclg);
+					}
+					sqfclg = 0.5 * sqfclg;
+
+
+					printf("Here place 2\n");
+					printf("ix[0] = %lf\n",ix[0]);
+					printf("ix[1] = %lf\n",ix[1]);
+					printf("ix[2] = %lf\n",ix[2]);
+					printf("ix[3] = %lf\n",ix[3]);
+					printf("ix[4] = %lf\n",ix[4]);
+					printf("ix[5] = %lf\n",ix[5]);
+					printf("ix[6] = %lf\n",ix[6]);
+					printf("ix[7] = %lf\n",ix[7]);
+					printf("ix[8] = %lf\n",ix[8]);
+	
+					printf("sqdclg = %lf\n",sqfclg);
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+					
+
+
+				}
+// if (IB != 0) _>>>>
+			}else{
+				Cleb = sgnfac;
+				return Cleb;
+
+			}
+// if (IC < IB) _>>>
+
+		}else{
+
+			IT = IC;	
+			IC = IB; 
+			IB = IT; 
+			IT = IF; 
+			IF = -IE;
+			IE = -IT;	
+			sgnfac = sqrt((2.*A[2]+1.0)/(2.*A[1]+1.0));
+			if((int)(im[0]-im[3])/2 %2 != 0){sgnfac = (-1.)*sgnfac;}
+			
+			Cleb = sgnfac;
+			return Cleb; 
+		}
+	
+// if(IA < IB) _>>
+	}else{
+
+		if(IA >= IC){
+			IT = IC;
+			IC = IB;
+			IB = IT;
+			IT = IF; 	
+			IF = -IE; 
+			IE = -IT; 			
+			
+			Cleb = sgnfac; 
+			return Cleb; 
+// if (IA < IC) as well as (IA < IB)
+		}else{
+			
+			IT =IA;
+			IA =IB;
+			IB =IT;
+			IT =ID;
+			ID =IE; 
+			IE =IT; 
+			
+			if((int)(K1/2) %2 != 0) {sgnfac = -1.;}
+//call to 135 again.		
+			if(IB != 0){
+				if(IE < 0){
+				//assuming the logic passes;
+					FC2 = IC +1;
+
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+					printf("Here place 3\n");
+					printf("ix[0] = %lf\n",ix[0]);
+					printf("ix[1] = %lf\n",ix[1]);
+					printf("ix[2] = %lf\n",ix[2]);
+					printf("ix[3] = %lf\n",ix[3]);
+					printf("ix[4] = %lf\n",ix[4]);
+					printf("ix[5] = %lf\n",ix[5]);
+					printf("ix[6] = %lf\n",ix[6]);
+					printf("ix[7] = %lf\n",ix[7]);
+					printf("ix[8] = %lf\n",ix[8]);
+	
+					printf("sqdclg = %lf\n",sqfclg);
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+				}else{
+	// if (IE >=0) _>>>
+					ID = -ID;
+					IE = -IE;
+					IF = -IF;
+					if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1.)*sgnfac;
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+					printf("Here place 4\n");
+					printf("ix[0] = %lf\n",ix[0]);
+					printf("ix[1] = %lf\n",ix[1]);
+					printf("ix[2] = %lf\n",ix[2]);
+					printf("ix[3] = %lf\n",ix[3]);
+					printf("ix[4] = %lf\n",ix[4]);
+					printf("ix[5] = %lf\n",ix[5]);
+					printf("ix[6] = %lf\n",ix[6]);
+					printf("ix[7] = %lf\n",ix[7]);
+					printf("ix[8] = %lf\n",ix[8]);
+	
+					printf("sqdclg = %lf\n",sqfclg);
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+					
+
+
+				}
+// if (IB != 0) _>>>>
+			}else{
+				Cleb = sgnfac;
+				return Cleb;
+
+			}
+	
+
+
+
+
+
+
+
+	
+
+		}
+
+
+
+
+
+	} 
+
+
+
+	return Cleb;
+}
+
+
+double CG(double J, double M, double j1, double m1, double j2, double m2){
+	//recall that j1,m1 + j2,m2 = J,M
+
+	if(M != m1 + m2) return 0;
+
+	double Jmin  =  abs(j1 - j2);
+	double Jmax  =  j1+j2;
+
+	if(J < Jmin || Jmax < J) return 0;
+	
+	double a0  = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0);
+	double A0 = sqrt(a0);
+
+	double a = factorial(J+M)   *factorial(J-M);
+	double a1= factorial(j1+m1) *factorial(j1-m1);
+	double a2= factorial(j2+m2) *factorial(j2-m2);
+	double A = sqrt( a  * a1 * a2);
+
+	int pmax = min( min(j1+j2-J,j1-m1),j2 + m2);
+	
+	double cg = 0.;
+
+	for( int p =0; p<=pmax;p++){
+
+		double p1 = factorial(j1+j2-J-p);
+		double p2 = factorial(j1-m1-p);
+		double p3 = factorial(j2+m2-p);
+		double p4 = factorial(J -j2 +m1 +p);
+		double p5 = factorial(J -j1 -m2 +p);
+		double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5);
+			
+		cg += t;
+	}
+	return A0*A*cg;
+}
+
+double CG_a(double A[6]){
+	//recall that j1,m1 + j2,m2 = J,M
+	//testing assingments until the output is the same. 
+	double j1 = A[0];
+	double j2 = A[1]; 
+	double J = A[2]; 
+	double m1 = A[3];
+	double m2 = A[4]; 
+	double M = A[5];
+ 
+	
+ 	
+	if(M != m1 + m2) return 0;
+
+	double Jmin  =  abs(j1 - j2);
+	double Jmax  =  j1+j2;
+
+	if(J < Jmin || Jmax < J) return 0;
+	
+	double a0  = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0);
+	double A0 = sqrt(a0);
+
+	double a = factorial(J+M)   *factorial(J-M);
+	double a1= factorial(j1+m1) *factorial(j1-m1);
+	double a2= factorial(j2+m2) *factorial(j2-m2);
+	double A1 = sqrt( a  * a1 * a2);
+
+	int pmax = min( min(j1+j2-J,j1-m1),j2 + m2);
+	
+	double cg = 0.;
+
+	for( int p =0; p<=pmax;p++){
+
+		double p1 = factorial(j1+j2-J-p);
+		double p2 = factorial(j1-m1-p);
+		double p3 = factorial(j2+m2-p);
+		double p4 = factorial(J -j2 +m1 +p);
+		double p5 = factorial(J -j1 -m2 +p);
+		double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5);
+			
+		cg += t;
+	}
+	return A0*A1*cg;
+}
+double ThreeJsym(double j1, double m1, double j2, double m2, double j3, double m3){
+	//[j1 j2 j3] = (-1)^(j1-j2-m3)/ sqrt(2*j3+1) * CG(j3, -m3, j1, m1, j2, m2)
+	//[m1 m2 m3]
+	
+	
+	
+	return pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j3,-m3,j1,m1,j2,m2);
+}
+
+double SixJsym(double j1, double j2, double j3, double j4, double j5, double j6){
+
+//--------------------------------------------------------------------------------//
+	// The six j symbol describes the coupling between j1 j2 and j3 to J - j1.
+	// essentially a triangle of angular momentum rules between these. 
+	// j1 = j1
+	// j2 = j2
+	// j3 = j1 + j2
+	// j4 = j3
+	// j5 = J = j1 + j2 + j3
+	// j6 = j2 + j3
+// ----------------------------------------------------------------------------- //
+// the following conditions check the triangle selection rules
+
+	if( j3 < abs(j1 - j2) || j1 + j2 < j3) return 0; 
+	if( j6 < abs(j2 - j4) || j2 + j4 < j6) return 0; 
+	if( j5 < abs(j1 - j6) || j1 + j6 < j5) return 0; 
+	if( j5 < abs(j3 - j4) || j3 + j4 < j5) return 0; 
+// now that they have been checked, we can go ahead and calculate sixJ.
+
+	double sixj = 0.0;
+	float m1 = -j1;
+	float m2 = -j2;
+	float m3 = -j3;
+	float m4 = -j4;
+	float m5 = -j5;
+	float m6 = -j6;
+
+	for(; m1 <= j1; m1 = m1 +1){
+		for(; m2 <= j2; m2 = m2 +1){
+			for(; m3 <= j3; m3 = m3 + 1){
+				for(; m4 <= j4; m4 = m4 +1){
+					for(; m5 <= j5; m5 = m5 + 1){
+						for(; m6 <= j6; m6 = m6 +1){
+	
+							double h = (j1 - m1) + (j2 - m2) + (j3 -m3) + (j4 - m4) + (j5 - m5) + (j6 - m6);
+							double b1 = ThreeJsym(j1, -m1, j2, -m2, j3, -m3);
+							double b2 = ThreeJsym(j1, m1, j5, -m5, j6, m6); 
+							double b3 = ThreeJsym(j4, m4, j2, m2,  j6, -m6);
+							double b4 = ThreeJsym(j4, -m4, j5, m5, j3, m3);					
+							double b = b1 * b2 * b3 * b4;
+
+							sixj += pow(-1,h)*b;
+	
+						}
+					}
+				}
+			}
+		}
+	}
+
+	return sixj;
+}
+
+
+int main(int argc, char ** argv){
+
+
+	double A[6] = {1,1,2,0,0,0};
+
+
+	double cg = CG2(A);
+	double cgr = CG(1,1,2,0,0,0);
+	printf("------\n");
+	printf("CG = %lf\n",cg);
+	printf("CGR = %lf\n",cgr);
+
+
+	return 0; 
+}
+
+
+
diff --git a/Cleb.h b/Cleb.h
new file mode 100644
index 0000000..4d4be1e
--- /dev/null
+++ b/Cleb.h
@@ -0,0 +1,653 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+
+//#include "Coeff.h"
+
+using namespace std;
+
+
+
+double CG2(double A[6]){
+
+	double ii[6];
+	double im[6];
+	double ix[9];
+	double nx[5];
+
+	for(int i = 0;i<6;i++){
+		ii[i] = 2.0*A[i];
+	}
+	double Cleb = 0.;
+	double IA = ii[0];
+	double IB = ii[1];
+	double IC = ii[2];
+	double ID = ii[3];
+	double IE = ii[4];
+	double IF = ii[5];
+
+	if(ID + IE - IF != 0){return 0;}
+	double K0 = IA + IB + IC;
+	if((int)K0 % 2 !=0){return 0;}
+
+	double K1 = IA + IB - IC;
+	double K2 = IC + IA - IB;
+	double K3 = IB + IC - IA;
+	double K4 = IA - abs(IB-IC);
+	double K5 = IB - abs(IC-IA);
+	double K6 = IC - abs(IA-IB);
+
+	double karr[6]  = {K1,K2,K3,K4,K5,K6};
+	double kmn=karr[0];
+	double kmx=karr[0];
+	for(int i=0;i<6;i++)
+	{
+		if(kmn>karr[i])
+		{
+			kmn=karr[i];
+		}
+		else if(kmx<karr[i])
+		{
+			kmx = karr[i];
+		}
+	}
+	double K7 = kmn; 
+
+	printf("K1 = %lf\n",K1);
+	printf("K2 = %lf\n",K2);
+	printf("K3 = %lf\n",K3);
+	printf("K4 = %lf\n",K4);
+	printf("K5 = %lf\n",K5);
+	printf("K6 = %lf\n",K6);
+	printf("K7 = %lf\n",K7);
+
+
+	if(K7 <  0) return 0; 
+
+	for(int i = 0;i<3;i++){
+		if((int)(ii[i]+ii[i+3])%2 !=0)return 0;
+		if(ii[i] < abs(ii[i+3]))return 0; 
+	}
+	double sgnfac =1.0;
+
+	for(int i = 0; i<6;i++){
+		im[i] = ii[i];
+	}
+	double FC2, IT = 0.;
+	//-------------------------------------------------------------------------
+	if(IA >= IB){
+		if(IC >= IB){
+			if(IB != 0){
+				if(IE < 0){
+				//assuming the logic passes;
+					FC2 = IC +1;
+
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+				}else{
+	// if (IE >=0) _>>>
+					ID = -ID;
+					IE = -IE;
+					IF = -IF;
+					if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1)*sgnfac;
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+					
+
+
+				}
+// if (IB != 0) _>>>>
+			}else{
+				Cleb = sgnfac;
+				return Cleb;
+
+			}
+// if (IC < IB) _>>>
+
+		}else{
+
+			IT = IC;	
+			IC = IB; 
+			IB = IT; 
+			IT = IF; 
+			IF = -IE;
+			IE = -IT;	
+			sgnfac = sqrt((2.*A[2]+1.0)/(2.*A[1]+1.0));
+			if((int)(im[0]-im[3])/2 %2 != 0){sgnfac = (-1.)*sgnfac;}
+			
+			Cleb = sgnfac;
+			return Cleb; 
+		}
+	
+// if(IA < IB) _>>
+	}else{
+
+		if(IA >= IC){
+			IT = IC;
+			IC = IB;
+			IB = IT;
+			IT = IF; 	
+			IF = -IE; 
+			IE = -IT; 			
+			
+			Cleb = sgnfac; 
+			return Cleb; 
+// if (IA < IC) as well as (IA < IB)
+		}else{
+			
+			IT =IA;
+			IA =IB;
+			IB =IT;
+			IT =ID;
+			ID =IE; 
+			IE =IT; 
+			
+			if((int)(K1/2) %2 != 0) {sgnfac = -1.;}
+//call to 135 again.		
+			if(IB != 0){
+				if(IE < 0){
+				//assuming the logic passes;
+					FC2 = IC +1;
+
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+				}else{
+	// if (IE >=0) _>>>
+					ID = -ID;
+					IE = -IE;
+					IF = -IF;
+					if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1.)*sgnfac;
+
+					IT  = K0/2 + 1;
+					ix[0] = IT-IC;
+					ix[1] = IT-IB;
+					ix[2] = IT-IA;
+					ix[3] = (IA+ID)/2 +1;
+					ix[4] = ix[3] - ID;
+					ix[5] = (IB+IE)/2 + 1;
+					ix[6] = ix[5] - IE;
+					ix[7] = (IC+IF)/2 +1;
+					ix[8] = ix[7] - IF;  
+					//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
+					double sqfclg = log(FC2) - lgamma(IT+1);
+					double IXI;
+					for(int i = 0; i<9;i++){
+						IXI = ix[i];
+						sqfclg = sqfclg + lgamma(IXI);
+					}
+					sqfclg = 0.5 * sqfclg;
+
+					double xarr[3] = {ix[0],ix[4],ix[5]};
+
+					double xmn=xarr[0];
+					double xmx=xarr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(xmn>xarr[i])
+						{
+							xmn=xarr[i];
+						}
+						else if(xmx<xarr[i])
+						{
+							xmx = xarr[i];
+						}
+					}
+					double NZMAX = xmn;
+
+					double NZ2 = (IB-IC-ID)/2;
+					double NZ3 = (IA-IC+IE)/2;
+
+					double carr[3] = {0, NZ2, NZ3};
+
+					double cmn=carr[0];
+					double cmx=carr[0];
+					for(int i=0;i<3;i++)
+					{
+						if(cmn>carr[i])
+						{
+							cmn=carr[i];
+						}
+						else if(cmx<carr[i])
+						{
+							cmx = carr[i];
+						}
+					}
+
+
+					double NZMIN = cmx+1;
+
+					if(NZMAX > NZMIN){return 0;}
+					double SS=0.;
+					double S1= pow((-1.),(NZMIN-1));
+					for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
+
+						double NZM1 = NZ -1;
+						nx[0] = ix[0]-NZM1;
+						nx[1] = ix[4]-NZM1;
+						nx[2] = ix[5]-NZM1;
+						nx[3] = NZ-NZ2;
+						nx[4] = NZ-NZ3;
+						double termlg = sqfclg - lgamma(NZ);
+						for(int i = 0; i<5; i++){
+							IXI = nx[i];
+							termlg= termlg -lgamma(IXI);
+							SS = SS + S1*exp(termlg);
+							S1 =-S1;
+							Cleb=sgnfac*SS;
+						}
+					} 
+					
+
+
+				}
+// if (IB != 0) _>>>>
+			}else{
+				Cleb = sgnfac;
+				return Cleb;
+
+			}
+	
+
+
+
+
+
+
+
+	
+
+		}
+
+
+
+
+
+	} 
+
+
+
+	return Cleb;
+}
+
+
+double CG(double J, double M, double j1, double m1, double j2, double m2){
+	//recall that j1,m1 + j2,m2 = J,M
+
+	if(M != m1 + m2) return 0;
+
+	double Jmin  =  abs(j1 - j2);
+	double Jmax  =  j1+j2;
+
+	if(J < Jmin || Jmax < J) return 0;
+	
+	double a0  = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0);
+	double A0 = sqrt(a0);
+
+	double a = factorial(J+M)   *factorial(J-M);
+	double a1= factorial(j1+m1) *factorial(j1-m1);
+	double a2= factorial(j2+m2) *factorial(j2-m2);
+	double A = sqrt( a  * a1 * a2);
+
+	int pmax = min( min(j1+j2-J,j1-m1),j2 + m2);
+	
+	double cg = 0.;
+
+	for( int p =0; p<=pmax;p++){
+
+		double p1 = factorial(j1+j2-J-p);
+		double p2 = factorial(j1-m1-p);
+		double p3 = factorial(j2+m2-p);
+		double p4 = factorial(J -j2 +m1 +p);
+		double p5 = factorial(J -j1 -m2 +p);
+		double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5);
+			
+		cg += t;
+	}
+	return A0*A*cg;
+}
+
+double CG_a(double A[6]){
+	//recall that j1,m1 + j2,m2 = J,M
+	//testing assingments until the output is the same. 
+	double j1 = A[0];
+	double j2 = A[1]; 
+	double J = A[2]; 
+	double m1 = A[3];
+	double m2 = A[4]; 
+	double M = A[5];
+ 
+	if(M != m1 + m2) return 0;
+
+	double Jmin  =  abs(j1 - j2);
+	double Jmax  =  j1+j2;
+
+	if(J < Jmin || Jmax < J) return 0;
+	
+	double a0  = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0);
+	double A0 = sqrt(a0);
+
+	double a = factorial(J+M)   *factorial(J-M);
+	double a1= factorial(j1+m1) *factorial(j1-m1);
+	double a2= factorial(j2+m2) *factorial(j2-m2);
+	double A1 = sqrt( a  * a1 * a2);
+
+	int pmax = min( min(j1+j2-J,j1-m1),j2 + m2);
+	
+	double cg = 0.;
+
+	for( int p =0; p<=pmax;p++){
+
+		double p1 = factorial(j1+j2-J-p);
+		double p2 = factorial(j1-m1-p);
+		double p3 = factorial(j2+m2-p);
+		double p4 = factorial(J -j2 +m1 +p);
+		double p5 = factorial(J -j1 -m2 +p);
+		double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5);
+			
+		cg += t;
+	}
+	return A0*A1*cg;
+}
+double ThreeJsym(double j1, double m1, double j2, double m2, double j3, double m3){
+	//[j1 j2 j3] = (-1)^(j1-j2-m3)/ sqrt(2*j3+1) * CG(j3, -m3, j1, m1, j2, m2)
+	//[m1 m2 m3]
+	
+	
+	
+	return pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j3,-m3,j1,m1,j2,m2);
+}
+
+double SixJsym(double j1, double j2, double j3, double j4, double j5, double j6){
+
+//--------------------------------------------------------------------------------//
+	// The six j symbol describes the coupling between j1 j2 and j3 to J - j1.
+	// essentially a triangle of angular momentum rules between these. 
+	// j1 = j1
+	// j2 = j2
+	// j3 = j1 + j2
+	// j4 = j3
+	// j5 = J = j1 + j2 + j3
+	// j6 = j2 + j3
+// ----------------------------------------------------------------------------- //
+// the following conditions check the triangle selection rules
+
+	if( j3 < abs(j1 - j2) || j1 + j2 < j3) return 0; 
+	if( j6 < abs(j2 - j4) || j2 + j4 < j6) return 0; 
+	if( j5 < abs(j1 - j6) || j1 + j6 < j5) return 0; 
+	if( j5 < abs(j3 - j4) || j3 + j4 < j5) return 0; 
+// now that they have been checked, we can go ahead and calculate sixJ.
+
+	double sixj = 0.0;
+	float m1 = -j1;
+	float m2 = -j2;
+	float m3 = -j3;
+	float m4 = -j4;
+	float m5 = -j5;
+	float m6 = -j6;
+
+	for(; m1 <= j1; m1 = m1 +1){
+		for(; m2 <= j2; m2 = m2 +1){
+			for(; m3 <= j3; m3 = m3 + 1){
+				for(; m4 <= j4; m4 = m4 +1){
+					for(; m5 <= j5; m5 = m5 + 1){
+						for(; m6 <= j6; m6 = m6 +1){
+	
+							double h = (j1 - m1) + (j2 - m2) + (j3 -m3) + (j4 - m4) + (j5 - m5) + (j6 - m6);
+							double b1 = ThreeJsym(j1, -m1, j2, -m2, j3, -m3);
+							double b2 = ThreeJsym(j1, m1, j5, -m5, j6, m6); 
+							double b3 = ThreeJsym(j4, m4, j2, m2,  j6, -m6);
+							double b4 = ThreeJsym(j4, -m4, j5, m5, j3, m3);					
+							double b = b1 * b2 * b3 * b4;
+
+							sixj += pow(-1,h)*b;
+	
+						}
+					}
+				}
+			}
+		}
+	}
+
+	return sixj;
+}
+
+
diff --git a/Coeff.h b/Coeff.h
new file mode 100644
index 0000000..de6ca97
--- /dev/null
+++ b/Coeff.h
@@ -0,0 +1,71 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+
+using namespace std;
+
+
+int param_run(int dt, int gt, int at, int st, int ft, int jt){
+	int param = 0; 
+	if( dt == 1 && gt == 1 && at == 1 && st == 1 && ft == 1 && jt == 1){
+		param =0; 
+	}else param = 1; 
+
+	return param;
+}
+
+
+
+int StrToInt(std::string const& s)
+{
+    std::istringstream iss(s);
+    int value;
+
+    if (!(iss >> value)) throw std::runtime_error("invalid int");
+
+    return value;
+}
+
+
+int factorial(int fact){
+	for(int i=1;i<=fact;i++){
+		fact=fact*i;
+	} 
+	return fact;
+}
+
+
+double FACTLOG(int num){
+
+	double faclog[170];
+	int RI = 0;
+
+	faclog[0] = 0.;
+	faclog[1] = 0.;
+	
+	for( int i = 3; i < 170; i++){
+		RI = i - 1; 
+		faclog[i] = log(RI) + faclog[i-1];
+
+	}
+	
+	double flog = faclog[num];
+
+
+	return flog;
+}
+
+
+
+
+
+
+
diff --git a/GUI.h b/GUI.h
new file mode 100644
index 0000000..dea6673
--- /dev/null
+++ b/GUI.h
@@ -0,0 +1,422 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+// VERSION : 1.0
+// Created with thanks to Daniel Folds Holt - University of Cambridge Cavendish Lab
+
+
+
+
+using namespace std;
+
+
+class HistoGUI{
+
+
+	public:
+
+		HistoGUI(){}
+
+		Display * disp;
+		Window    wind;
+		XEvent    evt;
+		int       screen;
+
+		XColor xcolour;
+		XColor xcolour_red;
+		XColor xcolour_veryred;
+		Colormap cmap;
+
+
+		std::vector<double> x;
+		std::vector<double> y;
+		double max_x;
+		double min_x;
+		double max_y;
+		double min_y;
+
+		double pos_x;
+		double pos_y;
+		double new_pos_x;
+		double new_pos_y;
+
+		unsigned int width;  
+		unsigned int height;  
+		double window_width;
+		double window_height;
+		double width_scale; 
+		double x_offset;
+		double height_scale;
+		double y_offset;
+
+		double old_xl, old_xh, old_yl, old_yh;
+		double old_mouse_x, old_mouse_y;
+
+		int Init();
+		int SetData(std::vector<double> a, std::vector<double>b);
+		int Loop();	
+		void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); }
+		int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win);
+		int DrawCrosshairs(int mouse_x, int mouse_y);
+		int Zoom(int mouse_x, int mouse_y);
+
+
+//		int Legendre_fit(std::vector<double> d1, std::vector<double> d2, double 
+
+};
+
+int HistoGUI::SetData(std::vector<double> a, std::vector<double>b){
+	for(int i=0; i < a.size(); i++){
+		x.push_back(a[i]);
+		y.push_back(b[i]);
+	}
+	return a.size();
+}
+
+int HistoGUI::Init(){
+
+	disp = XOpenDisplay(NULL);
+	if (disp == NULL) {
+		fprintf(stderr, "Cannot open display\n");
+		exit(1);
+	}
+
+	screen = DefaultScreen(disp);
+	wind   = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1,
+			BlackPixel(disp, screen), WhitePixel(disp, screen));
+	XGrabPointer(disp, wind, False, ButtonPressMask | ButtonReleaseMask | PointerMotionMask, GrabModeAsync,
+			GrabModeAsync, None, None, CurrentTime);
+	XSelectInput(disp, wind, ExposureMask | KeyPressMask | ButtonPressMask | ButtonReleaseMask | PointerMotionMask);
+	XMapWindow(disp, wind);
+
+	// colours
+	cmap = DefaultColormap(disp, screen);
+	xcolour.red = 32000; 
+	xcolour.green = 32000; 
+	xcolour.blue = 42000;
+	xcolour.flags = DoRed | DoGreen | DoBlue;
+	XAllocColor(disp, cmap, &xcolour);
+
+	xcolour_red.red   = 42000; 
+	xcolour_red.green = 32000; 
+	xcolour_red.blue  = 32000; 
+	xcolour_red.flags = DoRed | DoGreen | DoBlue;
+	XAllocColor(disp, cmap, &xcolour_red);
+
+	xcolour_veryred.red   = 65000; 
+	xcolour_veryred.green = 32000; 
+	xcolour_veryred.blue  = 32000; 
+	xcolour_veryred.flags = DoRed | DoGreen | DoBlue;
+	XAllocColor(disp, cmap, &xcolour_veryred);
+
+
+
+	return 1;
+
+}
+
+int HistoGUI::DrawCrosshairs(int mouse_x, int mouse_y){
+
+	int j1, j2;
+	unsigned int j3, j4;  
+	Window root_return;
+
+	XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4);
+	XSetForeground(disp,  DefaultGC(disp, screen), xcolour_red.pixel);
+	XDrawLine(disp, wind, DefaultGC(disp, screen), mouse_x, 0.0, mouse_x, height);
+	XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, mouse_y, width,  mouse_y);
+
+	pos_y = mouse_y * height_scale - y_offset;
+	pos_x = mouse_x * width_scale  - x_offset;
+
+	char coord[32];
+	sprintf(coord, "(%.2f, %.2f)", pos_x, pos_y);
+	XDrawString(disp, wind, DefaultGC(disp, screen), mouse_x + 5, mouse_y + 12, coord, strlen(coord));
+
+	return 1;
+}
+
+int HistoGUI::Zoom(int mouse_x, int mouse_y){
+
+	new_pos_x = mouse_x * width_scale - x_offset;
+	new_pos_y = mouse_y * height_scale  - y_offset;
+
+	pos_x = old_mouse_x * width_scale - x_offset;
+	pos_y = old_mouse_y * height_scale  - y_offset;
+
+
+	if(std::abs(mouse_x - old_mouse_x) < 20 or std::abs(mouse_y - old_mouse_y) < 20) return 1;
+
+
+	//printf("Zooming: %f - %f, %f - %f\n", mouse_x, old_mouse_x,mouse_y,old_mouse_y);
+
+	double x_low_win, y_low_win, x_hi_win, y_hi_win;
+	if(new_pos_x < pos_x){
+		x_low_win = new_pos_x;
+		x_hi_win  = pos_x;
+		//printf("new low; [%f, %f]\n",x_low_win, x_hi_win);
+	} else {
+		x_low_win = pos_x;
+		x_hi_win  = new_pos_x;
+		//printf("new high; [%f, %f]\n",x_low_win, x_hi_win);
+	}
+	if(new_pos_y > pos_y){
+		y_low_win = new_pos_y;
+		y_hi_win  = pos_y;
+	} else {
+		y_low_win = pos_y;
+		y_hi_win  = new_pos_y;
+	}
+
+
+	//printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win);
+	XClearWindow(disp, wind);
+	DrawData(x_low_win, y_low_win, x_hi_win, y_hi_win);
+	//double x_low_win, double y_low_win, double x_hi_win, double y_hi_win
+
+	return 1;
+
+}
+
+
+int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){
+
+
+	int j1, j2;
+	unsigned int j3, j4;  
+	Window root_return;
+
+	XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4);
+
+	//printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2);
+	//printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win);
+
+	double x_step;// = (width  * 0.8) / x.size();
+	double y_step;// = (height * 0.8) / x.size();
+
+
+	if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){	
+		// Audomatically decide data postition
+		max_x = x[0];
+		max_y = y[0];
+		min_x = x[0];
+		min_y = y[0];
+
+		for(int i=0; i<x.size(); i++){
+			if(x[i] > max_x) max_x = x[i];
+			if(x[i] < min_x) min_x = x[i];
+			if(y[i] > max_y) max_y = y[i];
+			if(y[i] < min_y) min_y = y[i];
+		}
+		//printf(" max_x = %f\n", max_x );
+		//printf(" max_y = %f\n", max_y );
+		//printf(" min_x = %f\n", min_x );
+		//printf(" min_y = %f\n", min_y );
+
+		width_scale = (max_x - min_x) / (0.8 * width);
+		x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x);
+
+		height_scale = -1. * (max_y - min_y) / (0.8 * height);
+		y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y);
+
+		//printf("width scale  = %f\n", width_scale);
+		//printf("x_offset     = %f\n", x_offset);
+		//printf("height scale = %f\n", height_scale);
+		//printf("y_offset     = %f\n", y_offset);
+		double axis_x  = (0. + x_offset) / width_scale;
+		double axis_y  = (0. + y_offset) / height_scale;
+
+		XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y);
+
+		char axis_val[4];
+		int w_step = width / 10;
+		for(int i=0; i < (int) width; i += w_step){
+			double x_val   = i * width_scale - x_offset;
+			sprintf(axis_val, "%.1f", x_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10, axis_val, strlen(axis_val));
+		}
+
+		int h_step = height / 10;
+		for(int i=0; i < (int) height; i += h_step){
+			double y_val = i * height_scale - y_offset;
+			sprintf(axis_val, "%.1f", y_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10, i, axis_val, strlen(axis_val));
+		}
+
+
+		XSetForeground(disp, DefaultGC(disp,screen), 0);	
+		double x_wid ;
+		double y_wid ;
+		double x_wid2;
+		double y_wid2;
+
+		for(int i=0; i < x.size() - 2; i++){
+			x_wid  = (x[i] + x_offset) / width_scale;
+			y_wid  = (y[i] + y_offset) / height_scale;
+			x_wid2 = (x[i + 1] + x_offset) / width_scale;
+			y_wid2 = (y[i + 1] + y_offset) / height_scale;
+			//printf("(%f, %f), (%f,%f)\n", x_wid,y_wid,x_wid2,y_wid2);
+			XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2);
+			XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4);
+		}
+		XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4);
+
+	} else {
+		//double x_low_win, double y_low_win, double x_hi_win, double y_hi_win
+
+		width_scale = (x_hi_win - x_low_win) / width;
+		x_offset = -1.0 * x_low_win;
+
+		height_scale = (y_hi_win - y_low_win) / height;
+		y_offset = -1.0 * y_low_win;
+
+		//printf("width scale  = %f\n", width_scale);
+		//printf("x_offset     = %f\n", x_offset);
+		//printf("height scale = %f\n", height_scale);
+		//printf("y_offset     = %f\n", y_offset);
+		double axis_x  = (0. + x_offset) / width_scale;
+		double axis_y  = (0. + y_offset) / height_scale;
+
+		XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y);
+
+		char axis_val[4];
+		int w_step = width / 10;
+		for(int i=0; i < (int) width; i += w_step){
+			double x_val   = i * width_scale - x_offset;
+			sprintf(axis_val, "%.1f", x_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10, axis_val, strlen(axis_val));
+		}
+
+		int h_step = height / 10;
+		for(int i=0; i < (int) height; i += h_step){
+			double y_val = i * height_scale - y_offset;
+			sprintf(axis_val, "%.1f", y_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10, i, axis_val, strlen(axis_val));
+		}
+
+
+		XSetForeground(disp, DefaultGC(disp,screen), 0);	
+		double x_wid ;
+		double y_wid ;
+		double x_wid2;
+		double y_wid2;
+
+		for(int i=0; i < x.size() - 2; i++){
+			x_wid  = (x[i] + x_offset) / width_scale;
+			y_wid  = (y[i] + y_offset) / height_scale;
+			x_wid2 = (x[i + 1] + x_offset) / width_scale;
+			y_wid2 = (y[i + 1] + y_offset) / height_scale;
+			//	printf("(%f, %f), (%f,%f)\n", x_wid,y_wid,x_wid2,y_wid2);
+			XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2);
+			XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4);
+		}
+		XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4);
+
+
+	}
+
+	old_xl = x_low_win;
+	old_yl = y_low_win;
+	old_xh = x_hi_win;
+	old_yh = y_hi_win;
+
+	return 1;
+}
+
+
+int HistoGUI::Loop(){
+
+	bool MousePressed  = false;
+	bool MousePressed2 = false;
+
+	while (1) {
+		XNextEvent(disp, &evt);
+		if (evt.type == Expose) {
+			//		XFillRectangle(disp, wind, DefaultGC(disp, screen), 20, 20, 10, 10);
+			//			XDrawString   (disp, wind, DefaultGC(disp, screen), 10, 50, msg, strlen(msg));
+			DrawData(-1,-1,-1,-1);
+		} else if (evt.type == ButtonPress){
+			/* store the mouse button coordinates in 'int' variables. */
+			/* also store the ID of the window on which the mouse was */
+			/* pressed.                                               */
+			int mouse_x = evt.xbutton.x;
+			int mouse_y = evt.xbutton.y;
+			old_mouse_x = mouse_x;
+			old_mouse_y = mouse_y;
+
+			/* check which mouse button was pressed, and act accordingly. */
+			if(evt.xbutton.button == Button1){
+				/* draw a pixel at the mouse position. */
+				XClearWindow(disp, wind);
+				DrawCrosshairs(mouse_x, mouse_y);
+				DrawData(old_xl, old_yl, old_xh, old_yh);
+				MousePressed = true;
+				printf("(%f, %f)\n", pos_x, pos_y);
+			} if(evt.xbutton.button == Button3){
+				MousePressed2 = true;
+			} if(evt.xbutton.button == Button2){
+				XClearWindow(disp, wind);
+				DrawData(old_xl, old_yl, old_xh, old_yh);
+			}
+		} else if (evt.type == ButtonRelease){
+			/* store the mouse button coordinates in 'int' variables. */
+			/* also store the ID of the window on which the mouse was */
+			/* pressed.                                               */
+			int mouse_x = evt.xbutton.x;
+			int mouse_y = evt.xbutton.y;
+			//printf("Released button: %i, %i\n", mouse_x, mouse_y);
+			//printf("Pressed button: %i, %i\n", mouse_x, mouse_y);
+			/* check which mouse button was pressed, and act accordingly. */
+			if(evt.xbutton.button == Button1){
+				/* draw a pixel at the mouse position. */
+				Zoom(mouse_x, mouse_y);
+				//DrawData();
+				MousePressed = false;
+			} if(evt.xbutton.button == Button3){
+				MousePressed2 = false;
+			}
+		} else if (evt.type == MotionNotify and MousePressed){
+
+			int mouse_x = evt.xmotion.x;
+			int mouse_y = evt.xmotion.y;
+			XClearWindow(disp, wind);
+
+			DrawCrosshairs(old_mouse_x, old_mouse_y);
+			DrawCrosshairs(mouse_x, mouse_y);
+			DrawData(old_xl, old_yl, old_xh, old_yh);
+
+		} else if (evt.type == MotionNotify and MousePressed2){
+
+			int mouse_x = evt.xmotion.x;
+			int mouse_y = evt.xmotion.y;
+
+			XSetForeground(disp, DefaultGC(disp,screen), xcolour_veryred.pixel);
+			XFillRectangle(disp, wind, DefaultGC(disp, screen), mouse_x -2, mouse_y -2, 4, 4);
+
+		} else if (evt.type == KeyPress){
+
+			printf("Key pressed = %x\n", evt.xkey.keycode);
+			if(evt.xkey.keycode == 0x41 or evt.xkey.keycode == 0x39){
+				XClearWindow(disp, wind);
+				DrawData(-1,-1,-1,-1);
+
+			} else {
+				break;	
+			}
+		}
+	}
+
+	return 1;
+}
diff --git a/QDK.h b/QDK.h
new file mode 100644
index 0000000..b915169
--- /dev/null
+++ b/QDK.h
@@ -0,0 +1,206 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+
+using namespace std;
+
+
+int QKN(double Energy, double radius, double distance, double thickness){
+	
+	
+	double Qkn = 0.; 
+
+	double E_mev = Energy/1000;
+	
+	double E_log = log(E_mev);
+
+	double EL1 = E_log;
+	double EL2 = pow(E_log,2);
+	double EL3 = EL1*EL2;
+	double EL4 = pow(EL2,2);
+	double EL5 = EL4*EL1;
+
+	double TT = -1.1907 -0.5372*EL1 - 0.0438*EL2 + 0.0218*EL3 + 0.0765*EL4 + 0.0095*EL5;
+
+	double Tau = exp(TT);
+
+	//calulating attenuation angles
+	
+	
+	double Z1 = radius / (distance + thickness);
+	double Z2 = radius / distance; 
+
+	double alpha = atan(Z1);
+	double gamma = atan(Z2);
+	double beta;
+
+
+	double BL = 0.;
+	double BU = alpha;
+	double A = 0.;
+	double delx1  = (BU-BL)/1000;
+	
+	double sum1,sum2,sum3 = 0.;
+	double sum4,sum5,sum6 = 0.;
+
+	double cosb,sinb,secb,c2,c4,fac1,ex1 = 0.;
+ 	double term1,term2,term3 = 0.;
+	double term4,term5,term6 = 0.;
+	int J=0;
+	
+	int loop_length = 500;
+
+	for(int i = 0; i<=loop_length; i++){
+		
+		if(i > 0 and i < loop_length){
+			J = i%2;	
+			//printf("\t\ti = %d\nJ=%d\n",i,J);
+			if(J==0){A=2.;
+			}else {A=4.;}
+			beta = BL+i+delx1;
+		}else{A=.1;beta = BL+i+delx1;}
+	
+		cosb = cos(beta);
+		sinb = sin(beta);
+		secb = 1.0/cosb;
+		c2 = pow(cosb,2);		
+		c4 = pow(cosb,4);
+		fac1 = -1 *Tau *thickness *secb;
+		ex1 = exp(fac1);
+
+		term1 = 0.5*(3*c2-1)*(1-ex1)*sinb*A*delx1;
+		term2 = 0.125*A*(35*c4-30*c2+3)*(1-ex1)*sinb*delx1;
+		term3 = A*(1-ex1)*sinb*delx1;
+
+		sum1 = sum1 +term1;
+		sum2 = sum2 +term2;
+		sum3 = sum3 +term3;
+
+	}
+		
+
+	double ans1 = sum1/3;
+	double ans2 = sum2/3;
+	double ans3 = sum3/3;	
+
+	double J2,B,beta2,cosb2,sinb2,secb2,cscb2,c22,c44,fac2,ex2 = 0.;
+	double LB=alpha;
+	double UB=gamma;
+	double delx2 = (UB-LB)/1000;
+
+	for(int i = 0; i<=loop_length; i++){
+		
+		if(i > 0 and i < loop_length){
+			J2 = i%2;	
+			//printf("\t\ti = %d\nJ=%d\n",i,J2);
+			if(J2==0){B=2.;
+			}else {B=4.;}
+			beta2 = LB+i+delx2;
+		}else{B=.1;beta2 = LB+i+delx2;}
+		
+		cosb2 = cos(beta2);
+		sinb2 = sin(beta2);
+		secb2 = 1.0/cosb2;
+		cscb2 = 1.0/sinb2;
+		c22 = pow(cosb2,2);		
+		c44 = pow(cosb2,4);
+		fac2 = -1 *Tau *(radius*cscb2 -distance*secb2);
+		ex2 = exp(fac2);
+/*
+		printf("cosb2---%lf\n",cosb2);
+		printf("sinb2---%lf\n",sinb2);
+		printf("secb2---%lf\n",secb2);
+		printf("cscb2---%lf\n",cscb2);
+		printf("c22---%lf\n",c22);
+		printf("c44---%lf\n",c44);
+		printf("fac2---%lf\n",fac2);
+		printf("ex2---%lf\n",ex2);
+		printf("-------------------------------------\n");
+*/
+		term4 = 0.5*(3*c22-1)*(1-ex2)*sinb2*B*delx2;
+		term5 = 0.125*B*(35*c44-30*c22+3)*(1-ex2)*sinb2*delx2;
+		term6 = B*(1-ex2)*sinb2*delx2;
+
+		sum4 = sum4 +term4;
+		sum5 = sum5 +term5;
+		sum6 = sum6 +term6;
+
+	}
+	double	ans4=sum4/3;
+	double	ans5=sum5/3;
+	double	ans6=sum6/3;
+
+/*
+	printf("ans1:%lf\n",ans1);
+	printf("ans2:%lf\n",ans2);
+	printf("ans3:%lf\n",ans3);
+	printf("ans4:%lf\n",ans4);
+	printf("ans5:%lf\n",ans5);
+	printf("ans6:%lf\n",ans6);
+*/	
+	double QD2  = (ans1+ans4)/(ans3+ans6);
+	double QD4  = (ans2+ans5)/(ans3+ans6);
+	printf("--------------\n");
+	printf("  QD2 = %lf\n",QD2);
+	printf("  QD4 = %lf\n",QD4);
+	printf("--------------\n");
+
+//Now output a file that contains R, D , T , gamma energy, attentuation coeff, q2 and q4
+
+	ofstream fileo;
+	fileo.open ("ad.txt");
+	fileo << "Radius = " << radius <<" [cm]\n";
+	fileo << "Distance = " << distance <<" [cm]\n";
+	fileo << "Thickness = " << thickness <<" [cm]\n";
+	fileo << "Atten.C = " << Tau <<" [cm^-1]\n";
+	fileo << "Gamma_E = " << Energy <<" [KeV]\n";
+	fileo << "QD2 = " << QD2 << "\n";
+	fileo << "QD4 = " << QD4 << "\n";
+	fileo.close();
+
+	return 1; 
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Racah.h b/Racah.h
new file mode 100644
index 0000000..2f7c85d
--- /dev/null
+++ b/Racah.h
@@ -0,0 +1,303 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+
+//#include "Coeff.h"
+
+using namespace std;
+
+
+
+double RACAH(double A[6]){
+	
+	double Racah = 0.;
+
+	double ii[6],ix[6],kk[8];
+	
+	double nn[4], nx[12],ny[8];
+
+	double K1,K2,K3,K4,K5,K6,K7,K8 = 0.;
+	
+	double NZMIN, NZMAX = .0;
+		
+	for(int i = 0;i<6;i++){
+		ii[i] = 2.0*A[i];
+	}
+	
+	double IA = ii[0];
+	double IB = ii[1];
+	double IC = ii[2];
+	double ID = ii[3];
+	double IE = ii[4];
+	double IF = ii[5];
+
+	K1 = IA + IB -IE;
+	K3 = IC + ID -IE;
+	K5 = IA + IC -IF;
+	K7 = IB + ID -IF;
+	K2 = IE - abs(IA-IB);
+	K4 = IE - abs(IC-ID);
+	K6 = IF - abs(IA-IC);
+	K8 = IF - abs(IB-ID);
+
+
+	double karr[8]  = {K1,K2,K3,K4,K5,K6,K7,K8};
+	double kmn=karr[0];
+	double kmx=karr[0];
+	for(int i=0;i<8;i++)
+	{
+		if(kmn>karr[i])
+		{
+			kmn=karr[i];
+		}
+		else if(kmx<karr[i])
+		{
+			kmx = karr[i];
+		}
+	}
+	double K0 = kmn; 
+	
+	if(K0 < 0) return 0;
+	
+	double nkk, n2kk,nkk2 = 0.;
+
+ 	kk[0] = K1;	
+ 	kk[1] = K2;	
+ 	kk[2] = K3;	
+ 	kk[3] = K4;	
+	kk[4] = K5;	
+ 	kk[5] = K6;	
+ 	kk[6] = K7;	
+ 	kk[7] = K8;	
+
+	for(int i = 0; i < 7; i = i + 2){
+		nkk = kk[i];
+		n2kk = kk[i]/2;
+		nkk2 = n2kk*2.;
+		if(nkk != nkk2){return 0;}
+
+	}
+	
+	ix[0] = IA;
+	ix[1] = ID;
+	ix[2] = IB;
+	ix[3] = IC;
+	ix[4] = IF; 
+	ix[5] = IE;
+	double Lmin = ix[0];
+	int Kmin = 1;
+	for( int i = 1; i < 6; i++){
+		if(ix[i] - Lmin < 0){
+			Lmin=ix[i];
+		}
+		if(ix[i] - Lmin == 0){
+			Kmin=i;
+		}
+		
+		if(ix[i] - Lmin > 0){
+			break;
+		}
+	}
+	
+	double sgn1 = 1.;
+	double KQ,IP,IQ,IT = 0.;
+	if(Kmin == 1){
+		KQ = 7 - Kmin;
+		
+		IA = ix[(int)KQ-1];
+		ID = ix[Kmin+4-1];
+		IE = ix[Kmin-1];
+		IF = ix[(int)KQ-4-1]; 
+
+	if(((int)(IA+ID-IE-IF)/2 %2) != 0) {sgn1 = -1.;}
+
+	}else if(Kmin == 2){
+	
+		KQ = 7 - Kmin;
+		
+		IA = ix[(int)KQ-1];
+		ID = ix[Kmin+4-1];
+		IE = ix[Kmin-1];
+		IF = ix[(int)KQ-4-1]; 
+
+		if(((int)(IA+ID-IE-IF)/2 %2) != 0) {sgn1 = -1.;}
+
+	}else if(Kmin == 3){
+		KQ = 9 - Kmin;
+		
+		IB = ix[(int)KQ-1];
+		IC = ix[Kmin+2-1];
+		IE = ix[Kmin-1];
+		IF = ix[(int)KQ -2-1];	
+		if(((int)(IB+IC-IE-IF)/2)%2 !=0) {sgn1 = -1.;}
+	}else if(Kmin == 4){
+		
+		KQ = 9 - Kmin;
+		
+		IB = ix[(int)KQ-1];
+		IC = ix[Kmin+2-1];
+		IE = ix[Kmin-1];
+		IF = ix[(int)KQ -2-1];	
+		if(((int)(IB+IC-IE-IF)/2)%2 !=0) {sgn1 = -1.;}
+	
+	}else if(Kmin == 5){
+
+		IE = ix[4];
+		IF = ix[5];
+		IB = ix[3];
+		IC = ix[2];
+
+	}else if(Kmin == 6){
+		
+	}
+	
+	IP = IA-IB;
+	IQ = IC-ID; 
+	
+	if(abs(IP) >= abs(IQ)){
+		if(IP >= 0){
+			IQ = IC-ID;
+		}else{
+
+		IT = IB;
+		IB = IA; 
+		IA = IT; 
+		IT = ID; 
+		ID = IC; 
+		IC = IT; 
+		IP = -IP;
+
+		}
+
+	}else{
+		IT = IC;
+		IC = IA; 
+		IA = IT; 
+		IT = ID;
+		ID = IB;
+		IB = IT;
+		IP = IQ;
+		if(IP >= 0){
+			IQ = IC - ID; 
+		}else{
+
+		IT = IB;
+		IB = IA; 
+		IA = IT; 
+		IT = ID; 
+		ID = IC; 
+		IC = IT; 
+		IP = -IP;
+
+		}
+	}
+
+	double sgn2 = 1.;	
+	
+	if((int)(IB+ID-IF)/2 %2 != 0) sgn2 = -1.;
+	if(IE >= 0){
+		double BI = IB; 
+		double DI = ID; 
+		Racah = sgn1 * sgn2 / sqrt((BI+1.)*(DI+1.));
+	}else{
+	
+	nn[0] = (IA+IB+IE)/2 +1;
+	nn[1] = (IC+ID+IE)/2 +1;
+	nn[2] = (IA+IC+IF)/2 +1;
+	nn[3] = (IB+ID+IF)/2 +1;
+	
+	nx[0] = nn[0] -IE;
+	nx[1] = nn[0] -IB;
+	nx[2] = nn[0] -IA;
+	nx[3] = nn[1] -IE;
+	nx[4] = nn[1] -ID;
+	nx[5] = nn[1] -IC;
+	nx[6] = nn[2] -IF;
+	nx[7] = nn[2] -IC;
+	nx[8] = nn[2] -IA;
+	nx[9] = nn[3] -IF;
+	nx[10] = nn[3] -ID;
+	nx[11] = nn[3] -IB;
+	
+	IP = (IA+IB+IC+ID+4)/2;
+	IQ = (IE+IF-IA-ID)/2;
+	IT = (IE+IF-IB-IC)/2;
+
+	double a = 1; 
+	double b = -IQ+1;
+	if(a > b) {
+
+		NZMIN = a; 
+	}else{
+		NZMIN = b;
+	
+	}
+	
+	double kzar[3]  = {nx[1],nx[4],nx[10]};
+	double kzmn=kzar[0];
+	double kzmx=kzar[0];
+	for(int i=0;i<3;i++)
+	{
+		if(kzmn>kzar[i])
+		{
+			kzmn=kzar[i];
+		}
+		else if(kzmx<kzar[i])
+		{
+			kzmx = kzar[i];
+		}
+	}
+	NZMAX = kzmn;
+ 
+	if(NZMAX < NZMIN){return 0;}
+	double sqlog = 0;
+	double I1 = 0.;
+	for(int i = 0; i <4; i++){
+		I1 = nn[i]+1;
+		sqlog = sqlog -lgamma(I1);
+	}
+	for(int i = 0; i <12; i++){
+		I1 = nx[i];
+		sqlog = sqlog - lgamma(I1);
+	}
+	sqlog = 0.5*sqlog;
+
+	double sss = 0.;
+	double NZM1 = 0.;
+	double sslog = 0.;
+
+	for(int nz = NZMIN; nz <= NZMAX; nz++){
+		NZM1 = nz - 1; 
+		ny[0] = IP - NZM1; 
+		ny[1] = nx[0] -NZM1;
+		ny[2] = nx[3] -NZM1;
+		ny[3] = nx[6] -NZM1;
+		ny[4] = nx[9]-NZM1;
+		ny[5] = nz;
+		ny[6] = IQ+nz;
+		ny[7] = IT+nz;
+
+		I1 = ny[0];
+		
+		sslog = sqlog + lgamma(I1);
+		for(int i = 1; i<8; i++){
+			I1 = ny[i];
+		}sslog = sslog - lgamma(I1);
+		sss = sss + (pow(-1.,NZM1)*exp(sslog));
+	}
+	Racah = sgn1*sss;	
+
+	}
+	return Racah;
+}
+	
+	
+
diff --git a/ad.txt b/ad.txt
new file mode 100644
index 0000000..de93b13
--- /dev/null
+++ b/ad.txt
@@ -0,0 +1,7 @@
+Radius = 3 [cm]
+Distance = 4 [cm]
+Thickness = 5 [cm]
+Atten.C = 0.294297 [cm^-1]
+Gamma_E = 1062 [KeV]
+QD2 = -0.499975
+QD4 = 0.374938