changing to use using complex class for alm

Author: Marion.Baumgartner

Holo_Multipol/CPPOut.dat

@@ -19,13 +19,13 @@ void Data::readData()
     dataF.open(infile);
     if(dataF.is_open())
     {
-        std::cout<<"reading file";
+        std::cout<<"--------------------------------------------------------------- \n reading file\n";
         for(int i=0;i<MAXANGLES;i++)
         {
-            double inNum0;
-            double inNum1;
-            double inNum2;
-            double inNum3;
+            double inNum0; // g(theta,phi)
+            double inNum1; // theta
+            double inNum2; //phi
+            double inNum3; //DOmega
             //std::cout<<"\n read data line "<<i;
             this->dataF>>inNum0;
             //std::cout<<" entry 0 "<<inNum0;
@@ -43,6 +43,7 @@ void Data::readData()
             {
                 thmax=inNum1;
             }
+            //convert angles from degres to radian
             messg[i][1]=deg_to_rad(inNum1);
             messg[i][2]=deg_to_rad(inNum2);
             messg[i][3]=deg_to_rad(inNum3);
@@ -56,12 +57,12 @@ void Data::readData()
         }
     }else
     {
-        std::cout<<"Can not open file! Will exit!";
+        std::cout<<"Can not open file! Will exit! \n ---------------------------------------------------------------";
         //Throw error here!!!
         return;
 
     }
-    std::cout<<"read file succesfully!";
+    std::cout<<"read file succesfully! \n ---------------------------------------------------------------";
     this->dataF.close();
 }
 

Holo_Multipol/Data.cpp

@@ -5,13 +5,13 @@
 Multipole::Multipole(std::string fileName, int lmax, int isym, double ekin): Data(fileName),LMAX(lmax), ISYM(isym), k(2*M_PI*sqrt(ekin/150))
 {
     assert(LMAX<Multipole::MAX_COEFF);
-    alm1.resize( LMAX/2+1, std::vector<double>());
-    alm2.resize( LMAX/2+1, std::vector<double>());
+    alm.resize( LMAX/2+1, std::vector<std::complex<double> >());
 }
 Multipole::~Multipole()
 {
-    delete alm1;
-    delete alm2;
+    //TODO enentiuelle leichen löschen
+    //delete alm1;
+    //delete alm2;
 }
 
 const int Multipole::getLMAX()
@@ -21,9 +21,8 @@ const int Multipole::getLMAX()
 
 void Multipole::multpl()
 {
-    std::cout<<"\n in multipl function"<<std::endl;
-    double dtheta=2*M_PI/180.0;
-    double bnorm=1;
+    std::cout<<"\n --------------------------------------------------------------- \n expanding multipole coefficients"<<std::endl;
+    double bnorm = 1;
     for(int l=0;l<=LMAX;l+=2)
     {
         int il=l/2;
@@ -34,32 +33,31 @@ void Multipole::multpl()
         //alm2[il].push_back(std::vector<double>(il));
         for(int m=0;m<=l;m+=ISYM)
         {
-            double rint1=0;
-            double rint2=0;
+            std::complex<double> int_res=0;
 
+            //go through all angles alm=int(g(theta,phi)*conj(Y_lm(thea,phi)) dOmega
             for(int i=0;i<MAXANGLES;i++)
             {
                 double phi=messg[i][2];
                 double theta=messg[i][1];
                 double gmessg=messg[i][0];
+                double dOmega=messg[i][3];
+
+                //std::cout<<"i, phi, theta, g "<<i<<" "<<phi<<","<<theta<<", "<<gmessg<<std::endl;
+                int_res+=gmessg*std::conj(boost::math::spherical_harmonic<double>(l, m, theta, phi))*dOmega;
 
-                std::cout<<"i, phi, theta, g "<<i<<" "<<phi<<","<<theta<<", "<<gmessg<<std::endl;
-                //norm=sqrt(4.0*M_PI/(2.0*l+1.0))
-                rint1+=vorz(m)*gmessg*boost::math::spherical_harmonic_r<double>(l, (-1)*m, theta, phi)*sin(theta)*messg[i][3];
 
-                rint2+=vorz(m)*gmessg*boost::math::spherical_harmonic_i<double>(l, (-1)*m, theta, phi)*sin(theta)*messg[i][3];
                 //std::cout<<"real, imag "<<rint1<<", "<<rint2<<"\n";
 
 
             }
-            if(l==0 && m==0 && rint1>0.001)
+            if(l==0 && m==0 && std::real(int_res)>0.001)
             {
-                    bnorm=rint1*dtheta;
+                bnorm=std::real(int_res);
                 std::cout<<"bnorm="<<bnorm<<std::endl;
             }
-            alm1[il].push_back(rint1*dtheta/bnorm);
-            alm2[il].push_back(rint2*dtheta/bnorm);
-            std::cout<<l<<" "<<m<<" "<<alm1[il][m]<<std::endl;//<<" "<<alm2[il][m]<<std::endl;
+            alm[il].push_back(int_res/ bnorm);
+            std::cout<<l<<" "<<m<<" "<<alm[il][m]<<std::endl;//<<" "<<alm2[il][m]<<std::endl;
         }
         //std::cout<<std::endl;
     }
@@ -69,42 +67,27 @@ std::cout<<"alm: (l,m,real,imag)\n";
 
 void Multipole::expans()
 {
-    std::cout<<"in expans function\n";
+    std::cout<<"\n --------------------------------------------------------------- \n expanding function according to calculated coefficients \n";
+    calc.clear();
     for(int i=0;i<MAXANGLES;i++)
     {
-        //std::cout<<"i "<<i<<std::endl;
+        std::cout<<"i "<<i<<std::endl;
         double theta=messg[i][1];
         double phi=messg[i][2];
-        double suml=0;
-        //std::cout<<"(theta, phi)=("<<theta<<", "<<phi<<")\n";
-        for (int l=0;l<LMAX;l+=2)
-        {
-            //std::cout<<"l is "<<l<<std::endl;
-            int il=l/2;
-            suml+=alm1[il][0]*boost::math::spherical_harmonic_r <double> (l,0,theta,phi); //m=0
-            //std::cout<<"sum"<<sum<<std::endl;
-            for(int m=1;m<=l;m++)
-            {
-                //std::cout<<"in seccond for loop!";
-                //calculate RE(sum_m A_lm*Y_lm)=sum_m RE(A_lm)*RE(Y_lm)-IM(A_lm)*IM(Y_lm)
-                double real=alm1[il][m]*boost::math::spherical_harmonic_r <double> (l,m,theta, phi); //imaginary part
-                double imag=alm2[il][m]*boost::math::spherical_harmonic_i <double> (l,m,theta, phi); //real part
-                suml+=2*(real-imag); //
-                //std::cout<<"sum, i"<<sum<<", "<<l<<std::endl;
-            }
-        }
+        double g_theta_phi=intencity(theta,phi);
         std::vector<double> temp;
-        std::cout<<suml<<" "<<180/M_PI*theta<<" "<<180/M_PI*phi<<std::endl;
-        temp.push_back(suml);
+        std::cout<<g_theta_phi<<" "<<180/M_PI*theta<<" "<<180/M_PI*phi<<std::endl;
+        temp.push_back(g_theta_phi);
         temp.push_back(theta);
         temp.push_back(phi);//TODO if time: make this call better -> works for a start
+
         calc.push_back(temp);
    }
-   std::cout<<"calc:(suml,theta,phi)\n";
+   std::cout<<"calc:(g(theat,phi),theta,phi)\n";
 }
 
 
-void Multipole::holorad(double alpha, double beta)
+/*void Multipole::holorad(double alpha, double beta)
 {
     double r=0;
     double kr=0;
@@ -243,25 +226,31 @@ void Multipole::smooth(double grid)
                 }
             }
         }
-    }*/
+    }*--->/
 }
 
 
 /*-------------------Private--------------------------------------*/
-double Multipole::innerSumm(int l, double alpha, double beta)
+double Multipole::intencity(double theta, double phi)
 {
-    double summ=0;
-    int il=l/2;
-    for(int m=ISYM;m<=l;m+=ISYM)
+    double summ0=0;
+    std::complex<double> summ (0,0);
+    for(int l=0; l<=LMAX; l+=2)
+    {
+        int m=0;
+        summ0+=std::real(alm[l/2][m]*boost::math::spherical_harmonic<double> (l,m,theta,phi));
+        m+=ISYM;
+        for(;m<=l;m+=ISYM)
         {
-            double real=alm1[il][m]*boost::math::spherical_harmonic_r <double> (l,m,alpha, beta); //imaginary part
-            double imag=alm2[il][m]*boost::math::spherical_harmonic_i <double> (l,m,alpha, beta);
-            summ+=2*(real-imag);
+            summ+=alm[l/2][m]*boost::math::spherical_harmonic<double> (l,m,theta,phi);
         }
-    return summ;
+    }
+
+
+    return summ0+2*std::real(summ);
 }
 
-void Multipole::printAlm()
+/*void Multipole::printAlm()
 {
     for(unsigned int i=0;i<alm1.max_size();i++)
     {
@@ -343,4 +332,4 @@ double Multipole::calcphi(double x, double y)
         }
     }
     return -1;
-}
+}*/

Holo_Multipol/Data.o

@@ -5,6 +5,8 @@
 #include <cmath>
 #include <vector>
 #include <assert.h>
+#include <complex>
+
 #include "Data.h"
 #include <boost/math/special_functions/legendre.hpp>
 #include <boost/math/special_functions/spherical_harmonic.hpp>
@@ -61,41 +63,41 @@ class Multipole:public Data //:public Data //this class now has acces to all the
     *The radial grid is in angstrom with 0.1\AA spacing
     *Function is adopted from the fortran Program by Juerg Osterwalde writen in 1993
     */
-    void holorad(double, double);
+    //void holorad(double, double);
 
     /**
     *Calculates a two dimensional image of the electron wave field near the photoemitter.
     *Equation (3) from A.Stucke et al. (1992) is used.
     *Ther radial grid is in angstroem with grid spacing in Angstroe
     *Function is adopted from the fortran Program by Juerg Osterwalde writen in 1993
     */
-    void doyzimage(double grid, int xyz);
+    //void doyzimage(double grid, int xyz);
 
     /**
     *Calculates a two-dimensional image of the electron wave field near the photoemitter.
     *Equation (3) from A.Stucke et al. (1992) is used.
     *the radial grid is in angstroem with grid angstroem spacing.
     *Function is adapted from the fortran program by Juerg Osterwalder written in 6.7.93
     */
-    void doxyzimage(double grid);
+    //void doxyzimage(double grid);
 
     /**
     *Takes a 2D image array and maps it into a 2d image array of integer numbers ranging from 0 to 255.
     *Function makes input in to 'image program more convenient for mac.
     *Calculates r * |img(xy)|**2 to produce images.
     *adapted from fortran program by Juerg Osterwalder, universite de fribourg, 12.7.93
     */
-    void scaleimage(double grid);
+    //void scaleimage(double grid);
 
     /**
     *Function performes a gausian smoothing of the image function
     */
-    void smooth(double grid);
+    //void smooth(double grid);
 
     /**
     *print the real coefficients that were calculated within this program
     */
-    void printAlm();
+    //void printAlm();
 
 //=======================================================================================//
     private:
@@ -104,8 +106,8 @@ class Multipole:public Data //:public Data //this class now has acces to all the
 
     const double k; //2*pi*sqrt(ekin/150)
 
-    std::vector<std::vector<double> > alm1; //real expansion coefficients
-    std::vector<std::vector<double> > alm2; //imaginary expansion coefficients
+    std::vector<std::vector<std::complex<double> > > alm; //expansion coefficients
+
     //To add elements to the two dimensional vector use alm.push_back(row)
     //where row is a vector of double type containing the data that belonges in to this row
     //std::vector<double> gCalc; take this vector from the base clas
@@ -114,23 +116,30 @@ class Multipole:public Data //:public Data //this class now has acces to all the
     *calculate (-1)^exp with exp in [1,2,3,4,...]
     *@param exp is a positiv or negative integer
     */
-    int vorz(int);
+    //int vorz(int);
 
 
-    inline double innerSumm(int l, double alpha, double beta);
+    /**
+    *calculate the intencity g(theta,phi) at a given point theta, phi using
+    *    g(theta,phi)=sum_{0<=l<=l_max}[A_l0*Y_l0(theta,phi))+2*Re(sum_{0<m<=l}A_lm*Y_lm(theta,phi))]
+
+    *@param theta    polar angle
+    *@param phi      azimutal angle
+    */
+    inline double intencity(double theta, double phi);
 
     /**
     *calculation of the proper polar angle within the yz-plane.
     *All polar angles are positive; negative y values should be accounted for by adding PI to the value of phi.
     *@return polar angle in radian (double precission)   ->  if the function returnes -1 there calculation was corrupded or the function is wrong
     */
-    double calcth(double y, double z);
+    //double calcth(double y, double z);
 
     /**
     *calculation of the proper azimutal angle within the yz-plane.
     *@return azimutal angle in radian (double precission)   ->  if the function returnes -1 there calculation was corrupded or the function is wrong
     */
-    double calcphi(double x, double y);
+    //double calcphi(double x, double y);
 
     static const int MAX_COEFF=100;//100 is the maximum amunt of coefficients that can be calculated
 };

Holo_Multipol/Multipole.cpp

@@ -18,6 +18,7 @@
 import griding2 as gr2
 
 
+
 class Grid:
     #Global variables
     def __init__(self, gfile="oldinp.itp", nfile="newinp.itp"):
@@ -256,7 +257,7 @@ def cPPOut(self):
         with open(self.outputfile, 'w\n') as outFile:
             #outFile.write("#g(theta, phi), theta, phi, dphi (values are in degrees)\n")
             for i in range(len(self.grid)):
-                value=str("%f %f %f %f" %(self.grid[i][3],self.grid[i][0],self.grid[i][1],self.grid[i][2]))
+                value=str("%f %f %f %f" %(self.grid[i][3],self.grid[i][0],self.grid[i][1],self.grid[i][5],))
                 #print value
                 outFile.write(value+"\n")
            # outFile.write("dtheta=%f"%self.dtheta)
@@ -322,12 +323,21 @@ def main():
     newGrid.fitNewToOldGrid()
     newGrid.writeGrid("temp.dat")
     #newGrid.plotGrid()
+    newGrid.cPPOut()
+    print"----------------------------------------------------"
+    print "writing output for C++ code to ",newGrid.outputfile
+    print"----------------------------------------------------"
 
-    calc = gr2.Calc(newGrid.grid,50)
+    calc = gr2.Calc(newGrid.grid,10)
 
     calc.multi()
     calc.expand()
     calc.writeData("calc.dat")
+    
+    
+    #print('\a')
+    #sys.stdout.write('\a')
+    #sys.stdout.flush()
     """s=-1
     while True:
         s = int(raw_input('Type 0 for Fortan output; 1 for c++ output '))