#include #include #include #include #include #include #include #include #include /* #include "TAxis.h" #include "TGraph.h" #include "TMultiGraph.h" #include "TCanvas.h" #include "TApplication.h" #include "TStyle.h" #include "TPad.h" #include "TROOT.h" #include "TColor.h" #include "TFrame.h" #include "TVirtualPad.h" #include "Math/SpecFunc.h" #include "THistPainter.h" */ using std::cout; using std::endl; //struct that stores all the input parameters struct inputParameters{ int n; int l; int m; int r_points; int aze_points; int pola_points; int r_max; }; double getX(double r ,double theta, double phi){ return r*cos(phi) * sin(theta); } double getY(double r ,double theta, double phi){ return r*sin(phi) * sin(theta); } double getZ(double r ,double theta, double phi){ return r*cos(theta); } int indDen(int i,int j , int k,int r_points, int pola_points){ //returns the right position in the 1D flat array that stores our 3d density //i = aze_points, j= pola_points, k = r_points return i * (r_points * pola_points) + j * (r_points) + k; } int factorial(int n) { if (n < 0) return 0; else return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; } std::complex PSI(double r, double theta, double phi, int n, int l,int m){ //have to cast the qunatum numbers to double double n_d = (double) n; double l_d = (double) l; double m_d = (double) m; //Hydrogen wavefunction therefore a_0 = 1 double coeff = sqrt(pow((2.0/n_d),3)) * factorial(n-l-1) / (2.0 * n_d*(double) factorial(n+l)); //Compute the generalized (associated) Laguerre polynomial of degree n and order k. double laguerre = std::assoc_laguerre(n-l-1,2*l+1,2.0*r/n_d); std::complex i1(0.0,1.0); std::complex sphHarm = std::sph_legendre(l,m,theta) * exp(i1 * m_d * phi); return coeff * exp(-r/n) * pow(2.0*r/n,l_d) * laguerre * sphHarm; } std::vector genDensity(double r_max, int r_points, int aze_points, int pola_points,int n, int l, int m){ //generates the density store it in a 3D array polar angle, azetumal angle, r_points int npoints = r_points * aze_points * pola_points; //store all elements in a single flat array to ensure continousness of memory std::vector density(npoints); for(int i = 0; i < aze_points; i++){ for(int j = 0; j < pola_points; j++){ for(int k = 0; k < r_points; k++){ double theta = M_PI / (double) aze_points * i; double phi = 2* M_PI / (double) pola_points * j; double r = r_max/ (double) r_points * k; std::complex PSI_value = PSI(r,theta,phi,n,l,m); //check phi <-> theta density[indDen(i,j,k,r_points,pola_points)] = real( conj(PSI_value) * PSI_value) ; } } } return density; } std::vector calcSurface(std::vector density, double r_max, int r_points, int aze_points, int pola_points){ //sort the array to get the the 90% biggest value //as this is done inplace we declare a new vector to call sort on so the //original vector does not get destroyed std::vector sort_density(density.begin(),density.end()); int max_elements_in_rays = r_points; int nRays = aze_points * pola_points; std::vector nmbr_in_ray(pola_points * aze_points); std::vector surface(pola_points * aze_points); sort(sort_density.begin(),sort_density.end()); int npoints = r_points * aze_points * pola_points; //the factor 0.9 is convention int ind_contour_value = (int) ceil(npoints * 0.9); double contour_value = sort_density[ind_contour_value]; //now we go through every "ray"(same azethumal and same polar angle) and look for the closest //value to contour value. We the add that element to a new vector which is the surface of the //to be plotted shape for(int i = 0; i < aze_points; i++){ for(int j = 0; j < pola_points; j++){ int p = 0; bool marker1 = 0; bool marker2 = 0; for(int k = 0; k < r_points; k++){ if(density[indDen(i,j,k,r_points,pola_points)] > contour_value) { marker1 = 1; } else{ marker1 = 0; } if(marker1 != marker2){ double r = r_max/ (double) r_points * k; double theta = M_PI / (double) aze_points*i; double phi = 2 * M_PI / (double) pola_points*j; //get kartesian coordinates for the poalar coordinates double x = getX(r,theta,phi); double y = getY(r,theta,phi); double z = getZ(r,theta,phi); cout << x <<" " << y << " " << z << "\n"; surface[i,j,p] = k; p++; } marker2 = marker1; } } } return surface; } int main(int argc, char *argv[]){ //The simulation parameters double r_max=10; int r_points=600; int pola_points=400; int aze_points=400; //orbitals that are plotted // n 1 2 3 4 4 // l 0 1 1 1 3 // m 0 0 0 0 3 //The orbital int n = 4; int l = 3; int m = 3; std::vector density = genDensity(r_max, r_points, aze_points,pola_points,n,l,m); std::vector surface = calcSurface(density, r_max, r_points, aze_points,pola_points); return 0; }