pderosa888/Angular_Distribution
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Add files via upload

Browse Source
This commit is contained in:
Peter DeRosa 2022-06-06 12:35:44 -04:00 committed by GitHub
parent d09ec29e9d
commit 9c61b7a85a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 3156 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
1062.csv Normal file
View File

@ -0,0 +1,4 @@
angle,Y,Yerr
45,129,10
90,110,10
135,129,10
1 angle Y Yerr
2 45 129 10
3 90 110 10
4 135 129 10

750
AD.cxx Normal file
View File

@ -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

740
Cleb.cxx Normal file
View File

@ -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;
}

653
Cleb.h Normal file
View File

@ -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;
}

71
Coeff.h Normal file
View File

@ -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;
}

422
GUI.h Normal file
View File

@ -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;
}

206
QDK.h Normal file
View File

@ -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;
}

303
Racah.h Normal file
View File

@ -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;
}

7
ad.txt Normal file
View File

@ -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