Please interpret valgrind output

[alexBB]

My environment is Ubuntu 16.06 and the code is .cpp

I posted a few days earlier here because of memory problems and the user Cubbi suggested that I use valgrind.

That was "Execution time failure" http://www.cplusplus.com/forum/unices/205986/

This is my compile statement:

g++ -std=c++14 -g sphere_perspective.cpp

and execution statement:

valgrind ./a.out

This is what I got:

==18099== Memcheck, a memory error detector
==18099== Copyright (C) 2002-2015, and GNU GPL'd, by Julian Seward et al.
==18099== Using Valgrind-3.11.0 and LibVEX; rerun with -h for copyright info
==18099== Command: ./a.out
==18099== 
*************SOME COPIOUS BUT IRRELEVANT OUTPUT*************
==18099== Mismatched free() / delete / delete []
==18099==    at 0x4C2EDEB: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x415D12: RealWork::legTransformNONrotated(int, std::complex<double>**, std::complex<double>**&, PhiLimits*) (sphere_perspective.cpp:2530)
==18099==    by 0x415029: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2295)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099==  Address 0x5ea2dd0 is 0 bytes inside a block of size 8 alloc'd
==18099==    at 0x4C2E80F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x41580C: RealWork::legTransformNONrotated(int, std::complex<double>**, std::complex<double>**&, PhiLimits*) (sphere_perspective.cpp:2497)
==18099==    by 0x415029: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2295)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099== 
jj =  0
==18099== Mismatched free() / delete / delete []
==18099==    at 0x4C2EDEB: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x4165AD: RealWork::restoreImageNONrotatedGaussLegendre(int, int, std::complex<double>**, PhiLimits*) (sphere_perspective.cpp:2591)
==18099==    by 0x415055: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2296)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099==  Address 0x6019510 is 0 bytes inside a block of size 8 alloc'd
==18099==    at 0x4C2E80F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x41614B: RealWork::restoreImageNONrotatedGaussLegendre(int, int, std::complex<double>**, PhiLimits*) (sphere_perspective.cpp:2572)
==18099==    by 0x415055: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2296)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099== 
jj =  1
jj =  2
qCoefficientsjj =  3
jj =  4
******************IRRELEVANT OUTPUT***************
jj =  19
==18099== Invalid write of size 8
==18099==    at 0x4152B8: RealWork::legendreNONrotatedAnalysis(int, PhiLimits*) (sphere_perspective.cpp:2330)
==18099==    by 0x415073: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2298)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099==  Address 0x7dcdf60 is 0 bytes after a block of size 256 alloc'd
==18099==    at 0x4C2E80F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x41520F: RealWork::legendreNONrotatedAnalysis(int, PhiLimits*) (sphere_perspective.cpp:2324)
==18099==    by 0x415073: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2298)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099== 
==18099== Mismatched free() / delete / delete []
==18099==    at 0x4C2EDEB: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x41561B: RealWork::legendreNONrotatedAnalysis(int, PhiLimits*) (sphere_perspective.cpp:2362)
==18099==    by 0x415073: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2298)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099==  Address 0x67e7030 is 0 bytes inside a block of size 8 alloc'd
==18099==    at 0x4C2E80F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x4153A3: RealWork::legendreNONrotatedAnalysis(int, PhiLimits*) (sphere_perspective.cpp:2339)
==18099==    by 0x415073: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2298)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099== 
==18099== Invalid read of size 8
==18099==    at 0x415467: RealWork::legendreNONrotatedAnalysis(int, PhiLimits*) (sphere_perspective.cpp:2348)
==18099==    by 0x415073: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2298)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099==  Address 0x5ac9ef8 is 0 bytes after a block of size 504 alloc'd
==18099==    at 0x4C2E80F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==18099==    by 0x40AB86: Localizations::fillPhiBoundsNONrotated(int, PhiLimits*&) (sphere_perspective.cpp:3571)
==18099==    by 0x414847: RealWork::integrateNONrotatedGaussLegendre(int) (sphere_perspective.cpp:2226)
==18099==    by 0x40AF82: main (sphere_perspective.cpp:4126)
==18099== 
==18099== 
==18099== HEAP SUMMARY:
==18099==     in use at exit: 5,017,760 bytes in 2,654 blocks
==18099==   total heap usage: 62,869,822 allocs, 62,867,168 frees, 6,380,959,633 bytes allocated
==18099== 
==18099== LEAK SUMMARY:
==18099==    definitely lost: 38,048 bytes in 8 blocks
==18099==    indirectly lost: 4,907,008 bytes in 2,645 blocks
==18099==      possibly lost: 0 bytes in 0 blocks
==18099==    still reachable: 72,704 bytes in 1 blocks
==18099==         suppressed: 0 bytes in 0 blocks
==18099== Rerun with --leak-check=full to see details of leaked memory
==18099== 
==18099== For counts of detected and suppressed errors, rerun with: -v
==18099== ERROR SUMMARY: 401665 errors from 5 contexts (suppressed: 0 from 0)
alex@alex-HP-ProBook-6555b:~/GFORTRAN-APPS/SPHERE-DESIGN/SPHERE-CPP$ g++ -std=c++14 -g sphere_perspective.cpp

Now, RealWork is a class. The others are routine names.

I would appreciate if someone will help me to analyze the valgrind output. What is that 18099?

Thanks you.

[ne555]

> What is that 18099?
the process id.

[alexBB]

OK, thanks, but how about the rest of it?

[TheIdeasMan]

You might get more clues from the compiler if you compile with this as a minimum:

g++ -std+c++14 -Wall -Wextra -pedantic-errors -g sphere_perspective.cpp

I use these ones as well:
http://www.cplusplus.com/forum/general/183731/#msg899203

Just because it compiles with appalling minimal options (IMO), doesn't mean it's right. Warnings are your friend :+)

Try not to use new and delete, find some other way. Mind you with that much code, doing so could be a real pain.

main (sphere_perspective.cpp:4126)

You have at least that many LOC in main() ? One should be able to split that up, one way or another. How many classes do you have?


Good Luck !!

[ne555]

> Mismatched free() / delete / delete []
> Address 0x5ea2dd0 is 0 bytes inside a block of size 8 alloc'd
delete what you new,
delete[] what you new[],
free() what you {m,c,re}alloc(),
don't mix them.

> Invalid read/write of size 8
> Address 0x7dcdf60 is 0 bytes after a block of size 256 alloc'd
out of bounds access: index goes from 0 to size-1.
array[size] is invalid.


> Try not to use new and delete, find some other way.
agreed, for example, you may use std::vector.

[TheIdeasMan]

So this is a Topic from 1 year ago:

http://www.cplusplus.com/forum/unices/183536/

Plus other related posts going back to Jan 2015 - actually most all of the OP's posts seem to be about this subject. I have read all of them. There seems to be a lot of experimentation and confusion. If one had a distinct lack of knowledge with C++ (at least in the beginning) why would one undertake such an involved project?

I understand you are converting a gfortran program to c++, there are some things that come to mind with that:

Procedural and OOP are completely different paradigms, so could require quite a bit design thought at the beginning. It could be a lot more involved than translating the code in place.

You have a lot of code in one file - that's a worry. Could you let us look at the code on a hosting site like Github?

I advised about using the warning options for the compiler about 1 year ago, and today I am doing it again. I wonder how much of the other advice given has been taken on board?

[alexBB]

Thank you all. Yes, it all started with FORTRAN 2003, I think, but I converted it to c++ more than a year ago. It is a complicated problem of integrating special functions on 2-Sphere. I've tried a few methods of integration, perhaps half a dozen and so far nothing works or rather half of a double integral works but the other part does not. I have a friend who is a professional c++ programmer, I gave the code to him and he could not crack it either. That was one of the reasons I converted the code to c++, he refused to look at FORTRAN. I can work on it only 3 days a week because I have a job and it is 4-10s. In the US it is a 4 day weekend now, a bonus.

People published similar code and claimed that it worked. Some code is from Germany. They are all professors at major universities. I cannot even compile any of them. When I write emails to them either nobody responds or as one German professor said: "You can play with this code but don't count on it." I get 215 compile errors when I put his code in Blend for Visual Studio under W-10. I get more errors in Ubuntu g++. The friend whom I mentioned spent an evening and could not compile it either.

I can put my code in github but I doubt you will enjoy it. It's still in a low stage of development with many, many, very many dead routines scattered around.

This code has a peculiar significance for me. I must have it working, I cannot explain now why but it is an absolute MUST. If this code works it will be only the beginning. Now it is just "to make it work," then I foresee optimization, adjusting parameters which are zillion to make it work super fast, etc, etc, etc. I have even a feeling I may need to go back to FORTRAN for speed.

I actually cracked the problem for now. This code compiled and I ran it. There were no run time errors I've compiled and run it numerous times before but I wrote a new routine with a new method of integration and forgot basics. Now I will keep it in mind and be very careful.

I am not answering all points of criticism but please explain to me why I should use malloc instead of new (...) and free(...). I never used delete. I don't know how valgrind come up with that. I use new and free.

I tried to use vectors long time ago, I think @TheIdeasMan suggested them, although I am not sure. For some reason I could never make them working.

I AM VERY GRATEFUL TO EVERYONE, especially @TheIdeasMan.

Thanks, - Alex

[ne555]

> I never used delete. I don't know how valgrind come up with that. I use new and free.
Again:
if you allocate with new, then you must release with delete
if you allocate with new[], then you must release with delete[]
if you allocate with {m,c,re}alloc(), then you must release with free()
don't mix them.


> but please explain to me why I should use malloc instead of new (...) and free(...)
sorry, I was not precise. It wasn't a criticism against new, but against manual memory management. Avoid `{m,c,re}alloc()', `new' and `new[]'.


> I can put my code in github but I doubt you will enjoy it.
¿do you have any good reason for not wanting your code to be public?
you could receive better help if people have access to your code.

[alexBB]

Thank you ne555. I will give it a thought. There is too much trash in the code, abandoned routines which did not work or I found a better solution. But I will give it a thought. I applied for a patent with USPTO but it was all math and various diagrams. This code is to implement it. I will discuss it with some people I trust.

Now the new/delete combination is clear. I will change it. Thank you very much.

Thanks, - Alex

[TheIdeasMan]

There is too much trash in the code, abandoned routines which did not work or I found a better solution.

Use version control - then you can revert to old versions, create branches and the like. git seems to be a popular choice. Make sure you stage and commit all you changes, to avoid losing stuff. Lots of IDE's have a GUI interface to it, making it easy to use.

It sounds as though you are not using an IDE on Ubuntu - it will make you life a lot easier if you do.

Now the new/delete combination is clear. I will change it. Thank you very much.

Not sure if you yet understand - ideally get rid of all forms of new / delete. Instead use the STL. I know you mentioned that you tried to use std::vector but couldn't get it to work - can you show us what you tried?

I applied for a patent with USPTO but it was all math and various diagrams.

Are you sure no one has done this already?

[TheIdeasMan]

Can you show us the fortran code? And some documentation about the problem itself?

I was thinking that we could look at it from a design point of view.

[alexBB]

I cannot do anything today, possibly I can do it next weekend since now it is late and tomorrow I will be at work. The problem conceptually appears to be broken in stages. At this stage which is not the first one I need to expand an image on 2-sphere (a 2-D function) into a series of spherical harmonics. In order to make sure everything is correct I want to reconstruct that original function. I need to do it only once, then I will move to other stages. The spherical harmonics must be strictly normalized and that was the first stage to generate them. There are many problems in doing so since some algorithms are unstable. The "expansion" is essentially numerical integration. The problem he is that integration in one dimension is simple but with two dimensions the complexity simply explodes. Everything becomes much more difficult by perhaps a factor of 10.

>Are you sure no one has done this already?

Yes, someone's done this already. It is me. I got a first patent four years ago.

[TheIdeasMan]

Right, so it is a given that there is some serious maths to be encountered. I probably won't be able to help very much there - it's been too long since I did any normal math, never mind anything complicated :+D

But having read all of your posts in the last 2 years, it seems that you have had a series of problems with the c++ side of things - it is in that area some help might be offered. One of the problems I have come across before is that someone has reasonable amount of code (in your case at least 4,000 LOC) and this code has been worked on for a long time, but there are major problems with it (from a c++ or OOP point of view), so it becomes a big job to fix it - even a complete re-write. Hence a great deal of reluctance to change.

There are many problems in doing so since some algorithms are unstable.

Do you mean mathematically, or did you have a problem with the code itself?

Just on the Legendre Polynomials, apparently C++17 is going to have built in functions to do those. There is reference material available now, but I don't think that it is part of g++ yet. Maybe in 2 or 3 months time. There is sph_legendre near the bottom :+)
http://en.cppreference.com/w/cpp/numeric/special_math

Apparently there is an implementation in boost already:
http://www.boost.org/doc/libs/1_63_0/libs/math/doc/html/math_toolkit/sf_poly/sph_harm.html

Hope this helps :+)

[TheIdeasMan]

Actually, it does exist on my system (I have g++ 6.3.1) - it's in tr1, as it says in cppreference :+D

gcc version 6.3.1 20161221 (Red Hat 6.3.1-1) (GCC)

#include "tr1/cmath"
#include <iostream>

int main()
{
    // spot check for l=3, m=0
    double x = 1.2345;
    std::cout << "Y_3^0(" << x << ") = " << std::tr1::sph_legendre(3, 0, x) << '\n';

    // exact solution
    double pi = std::acos(-1);
    std::cout << "exact solution = "
              << 0.25*std::sqrt(7/pi)*(5*std::pow(std::cos(x),3)-3*std::cos(x))
              << '\n';
}

[TheIdeasMan]

Yes, someone's done this already. It is me. I got a first patent four years ago.

Can you tell us what the patent number is, or provide a link to it?

At this stage which is not the first one I need to expand an image on 2-sphere ...

Consider what software like Google Earth does: I don't know exactly, but it would have to project imagery to either a sphere or ellipsoid.

[lastchance]

If it's any use, the code below should compute the expansion in spherical harmonics on the surface of the unit sphere. For more accurate recovery of the original function you will need more terms in the series (increase LMAX), but it will take longer.

The numerical integration is very basic (2-d mid-ordinate rule), but it isn't the accuracy-determining part.
(EDITED SEVERAL TIMES)

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <cmath>
#include <complex>
using namespace std;

const int LMAX = 50;                             // max l in series expansion
complex<double>coeff [1+LMAX][1+2*LMAX];         // coefficients in the expansion: coeff[l][l+m] is multiplier of Y_lm (** NOTE**)

const double PI = 3.14159265358979323846;

//======== Function prototypes

double testFunction( double theta, double phi );                               // test function
void computeCoefficients( int ntheta, int nphi );                              // computes coefficients in expansion
complex<double> surfaceIntegral( int ntheta, int nphi, int l, int m );         // does surface integral over sphere
complex<double> seriesSum( double theta, double phi );                         // sum of the series
void outputCoefficients();                                                     // outputs coefficients in expansion
void checkSeries();                                                            // checks series sum against original
double doubleFactorial( int n );                                               // double factorial n!!
double ZLM( int l, int m, double x );                                          // (Renormalised) associated Legendre
double associatedLegendre( int l, int m, double x );                           // Associated Legendre polynomial
complex<double> sphericalHarmonic( int l, int m, double theta, double phi );   // Spherical harmonic

//==================

int main()
{
   const int NTHETA = 100;               // number of angle intervals to do the surface integration
   const int NPHI   = 200;

   cout << "Computing series expansion ...\n\n";
   computeCoefficients( NTHETA, NPHI ); 
// outputCoefficients();

   checkSeries();                         // check against original function
}

//==================

double testFunction( double theta, double phi )       
{
   return cos( theta ) * sin( phi );           
}

//==================

void computeCoefficients( int ntheta, int nphi )     // c_lm = Integral f Ylm* dA;    NOTE: c_lm is held in coeff[l][l+m]
{
   for ( int l = 0; l <= LMAX; l++ )
   {
      for ( int m = -l; m <= l; m++ ) coeff[l][l+m] = surfaceIntegral( ntheta, nphi, l, m );
   }
}

//==================

complex<double> surfaceIntegral( int ntheta, int nphi, int l, int m )
{
   double theta;                          // colatitude
   double phi;                            // azimuth
   double dtheta = PI / ntheta;           // angle increment
   double dphi = 2.0 * PI / nphi;      
   double dA;                             // element of area

   complex<double> *Ylm0 = new complex<double>[ntheta];    // precompute theta part
   for ( int i = 1; i <= ntheta; i++ )     
   {                                                                 
      theta = ( i - 0.5 ) * dtheta;
      Ylm0[i-1] = sphericalHarmonic( l, m, theta, 0.0 );
   }

   complex<double> *phiPart = new complex<double>[nphi];   // precompute phi part
   phiPart[0] = polar( 1.0, -m * 0.5 * dphi );
   complex<double> multiplier = polar( 1.0, -m * dphi );
   for ( int j = 1; j < nphi; j++ ) phiPart[j] = multiplier * phiPart[j-1];

   complex<double> integral = 0.0;
   for ( int i = 1; i <= ntheta; i++ )           // colatitude loop
   {                                                                 
      theta = ( i - 0.5 ) * dtheta;
      dA = sin( theta ) * dphi * dtheta;
      for ( int j = 1; j <= nphi; j++ )          // azimuth loop
      {
         phi = ( j - 0.5 ) * dphi;
         integral += testFunction( theta, phi ) * Ylm0[i-1] * phiPart[j-1] * dA;
      }
   }

   delete [] Ylm0;
   delete [] phiPart;
   return integral;
}

//==================

complex<double> seriesSum( double theta, double phi )
{
   complex<double> value = 0.0;
   for ( int l = 0; l <= LMAX; l++ )
   {
      for ( int m = -l; m <= l; m++ ) value += coeff[l][l+m] * sphericalHarmonic( l, m, theta, phi );
   }
   return value;
}

//==================

void outputCoefficients()
{
   for ( int l = 0; l <= LMAX; l++ )
   {
      for ( int m = -l; m <= l; m++ )
      {
         cout << "l=" << setw(4) << l << "  m=" << setw(4) << m
              << "  c_lm=" << showpos << scientific << setw(16) << setprecision(6) 
              << coeff[l][l+m].real() << "  " << coeff[l][l+m].imag() << endl;
      }
   }
}

//==================

void checkSeries()
{
   double theta, phi;

   char ans = 'y';
   while( ans =='y' || ans == 'Y' )
   {
      cout << "Input colatitude (0 <= theta <= PI) and azimuth (0 <= phi <= 2 PI) in RADIANS: ";
      cin >> theta >> phi;
      cout << "Original function = " << testFunction( theta, phi )     << endl;
      cout << "Sum of series     = " << seriesSum( theta, phi ).real() << endl;
      cout << endl;
      cout << "Any more? (y/n): ";
      cin >> ans;
   }
}

//================== 

double doubleFactorial( int n )
{
   double f = 1.0;
   while ( n > 1 )
   {
      f *= n;
      n -= 2;
   }
   return f;
}

//================== 

double ZLM( int l, int m, double x )
{
   if ( l < 0 || m > l ) return 0.0;

   // Special cases
   if      ( l == 0 ) return 1.0;                          // Z_00
   else if ( l == 1 ) 
   {
      if ( m == 0 ) return x;                              // Z_10
      else          return -sqrt( 1.0 - x * x );           // Z_11
   }

   // Everything else by recursion
   double Zlm, Zl1m, Zl2m;
   Zl2m = ZLM( 0, m, x );                                  // recursive function call for one of the above
   Zl1m = ZLM( 1, m, x );
   for ( int i = 2; i <= l; i++ )
   {
      if ( i == m )
      {
         Zlm = pow( 1 - x * x, 0.5 * m );
         if ( m %2 != 0 ) Zlm = -Zlm;
      }
      else
      {
         Zlm = 0.0;
         if ( m <= i - 1 ) Zlm += Zl1m * x * ( 2 * i - 1 );      // recursion
         if ( m <= i - 2 ) Zlm -= Zl2m * ( i + m - 1 );
         if ( m < i ) Zlm /= ( i - m );
      }

      Zl2m = Zl1m;
      Zl1m = Zlm;
   }
   return Zlm;
}

//================== 

double associatedLegendre( int l, int m, double x )
{
   if ( l < 0 || abs( m ) > l ) return 0.0;

   int ma = abs( m );
   double Plm = ZLM( l, ma, x );
   if ( m >= 0 ) 
   {
      Plm = Plm * doubleFactorial( 2 * m - 1 );
   }
   else
   {
      for ( int i = 1; i <= ma; i++ ) Plm = Plm * ( 2.0 * i - 1.0 ) / ( ( l - ma + i ) * ( l + i ) );
      if ( m %2 != 0 ) Plm = -Plm;
   }
   return Plm;
}

//================== 

complex<double> sphericalHarmonic( int l, int m, double theta, double phi )
{
   int ma = abs( m );
   if ( l < 0 || ma > l ) return 0.0;

   double norm = ( 2.0 * l + 1.0 ) / ( 4.0 * PI );
   for ( int i = 1; i <= ma; i++ )
   {
      norm = norm * ( ( 2.0 * i - 1.0 ) / ( l - ma + i ) ) * ( ( 2.0 * i - 1.0 ) / ( l + i ) );
   }
   norm = sqrt( norm );

   complex<double> Ylm = norm * ZLM( l, ma, cos( theta ) ) * polar( 1.0, m * phi );
   if ( m < 0 && m % 2 != 0 ) Ylm = -Ylm;

   return Ylm;
}

//==================  

Input colatitude (0 <= theta <= PI) and azimuth (0 <= phi <= 2 PI) in RADIANS:  0.5  0.7
Original function = 0.565354
Sum of series     = 0.561265

Any more? (y/n):  y
Input colatitude (0 <= theta <= PI) and azimuth (0 <= phi <= 2 PI) in RADIANS:  2.3  5.8
Original function = 0.309553
Sum of series     = 0.309358

Any more? (y/n):  n

[alexBB]

I much appreciate all contributions especially by @lastchance and @TheIdeasMan. Computation of normalized Associated Legendre Polynomials (ALP) has taken some of my time for sure but now it seems to be behind. I've gone through a number of stages. First you write routines that compute the corresponding values for lower indexes l and m. The corresponding formulas are all over the Internet, I read them in Wikipedia. They have to be normalized and with the lower values of indexes it is not hard. Elementary multiplications do the job. There is a book written by a Russian Belousov which has been translated to English. It is a book of tables of ALP which are properly normalized with indexes l <=56.

https://www.amazon.com/Tables-Normalized-Associated-Legendre-Polynomials-ebook/dp/B01DJDETCO/ref=sr_1_7?ie=UTF8&qid=1485039301&sr=8-7&keywords=belousov

So, you take recurrence formulas and go up the ladder. The recurrence formulas first must be changed for them to contain only normalized Polynomials. Those formulas are also available in Wikipedia but they are not normalized. Various sources have normalized formulas but their derivation is elementary. Belousov used normalized formulas but I am using different ones.

While you compute in the area of l < 56 you can compare them with the Belousov's table, higher than that you are on your own. There are ways to improve the confidence. You first compute Legendre Polynomials (not associated) with m=0 and use recurrence formulas to go to the right increasing index m step by step. Then you compute sectorial harmonics or rather the corresponding polynomials (m=l) and using recurrence formulas move from right to left and when you met the previous value you just computed you see if they coincide.

Some sources claim to have computational methods that compute ALPs up to a few thousand (l=2800 for instance). I kind of doubt it because the values of such function around the North Pole must be pretty much zero.

In short it is a non-trivial exercise altogether and requires a lot of time and patience. As I said I think I am past this stage. I now operate in the area of l<=30 but may need up to l <= 600 soon.

Also if you use the method I employed (normalizing recurrence formulas first), you avoid dealing with factorials and double factorials which is a pain.

@lastchance, thank you. It seems you compute for l<=50. Very cute. I have to look into it.

Anyway, again I appreciate it all. Patent number? My name is over there, and I think my address also. I will wait with it, sorry.

It is a fascinating area I'll tell you. All what we discussed so far is a small part of what I do.

Thanks, - A.

[lastchance]

Hello @alexBB

In my code you can increase LMAX to much larger values for more terms, but it will take commensurately longer and I wanted users not to be left wondering if their code was actually still doing anything.

On the way these functions were calculated you are correct that normalisation is a challenge. I originally calculated the actual associated Legendre polynomials, but it broke down somewhere between LMAX = 40 and 50 because the coefficients become too big. If you delve into my code you will find (see functions ZLM() and associatedLegendre()) that I actually work with a renormalised function - I've called it ZLM - that is actually Plm divided by the (2m-1)!! double factorial for positive m (with a slightly different form for negative m). In fact I never actually use the full associated Legendre functions to calculate the spherical harmonics - they are only left here so that I could check against tabulated solutions. In the spherical harmonic itself I have written and re-ordered the normalisation factor to compensate for the (2m-1)!! factor and ensure boundedness at all times.

As you note I have coded the functions by two-point recursion, with all formulae being lifted from the Wolfram sites:
http://mathworld.wolfram.com/SphericalHarmonic.html
http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html
(also cross-checked with Wikipedia). I checked my final functions with the chart calculators at
http://keisan.casio.com/exec/system/1180573407

For the numerical integration I used a very simple mid-ordinate rule on a regular theta-phi grid (which puts an excessive number of points near N and S poles). Your original error message implied that you were using Gauss-Legendre quadrature. I tried this, but it didn't change my results at all (I had a very smooth test function), so I left that out.

As far as your original memory problem is concerned I only actually use new and delete [] in the numerical integration part. I think they balance correctly here, but I'm sure that the professional programmers on this forum will inform me if they don't!

Thanks for the interesting post. I covered these functions years ago as an undergraduate student, but we only dealt with mathematical symbols and properties there - it was quite enlightening to see what a challenge they are to actually compute. There's lots of applications for them in many areas of physics, but I hadn't noted any in image processing before - I wish you luck with your endeavour.

[alexBB]

@lastchance, thank you. Your post is very interesting. I haven't used wolfram. Casio website computes Associated Legendre polynomials but they are not normalized. It has been helpful in the past. There are actually many publications which use or attempt to use Spherical Harmonics in visual imaging but they all do it differently, so I hope I am a trailblazer, sort of. Otherwise it is CMB research, geodesy too. I think it is geodesy where they want to compute very high frequency polynomials. Anyhow, I would like to get in touch with you, just for a virtual handshake at first. My email is alex bb alex @ live.com Remove spaces to use it. It is my junk email account, so even if someone harvested it, I don't care. The stuff I get in my inbox in that account is mind boggling anyway. I've been very busy and so far haven't had a chance to delve in your code but I will, next weekend for sure. Thank you very much for the code again. - Alex

P.S. For numerical integration I've used Simpson, Newton, Gauss-Legendre as you noticed. I still have massive problems with that all. I'll tell you later.

[TheIdeasMan]

Hi,

Here is a link to geodetic equations as used in Australia:

http://www.icsm.gov.au/gda/gda-v_2.4.pdf

Note the Redfearn formulae, used to convert between Lat & Long to Grid (plane) co-ordinates. Also the formulae to convert to/from Earth Centred Earth Fixed (ECEF) cartesian co-ordinates and Lat/Long.

The equations are relatively simple compared to what it sounds like you are doing. And they are all Ellipsoidal based, not Spherical.

Hope it helps :+)