sphere.h

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// spectral toolkit 
// copyright (c) 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//

#ifndef __sphere__
#define __sphere__

#include "spectral.h"
#include "rfft.h"
#include "fiber.h"
#include "group.h"

namespace spectral
{	

  /// Scalar multiple instance spherical harmonic transform.  Computes the two-dimensional 
  /// spherical harmonic transform of a real scalar field on the sphere. 
  /// The grid is defined by the parameter N, which is the number of latitude
  /// pointer from the equator to the pole.  The total number of latitude
  /// points from pole to pole is 2N, with the number of longitude points set to 4N.
  /// See Swartztrauber, P.N., The Vector Harmonic Transform Method for 
  /// Solving Partial Differential Equations in Spherical Geometry,
  /// <I>Monthly Weather Review</I>, <B>121</B> (1993), 3415-3437.
  ///
  /// The following example computes a projection of an arbitrary 3-d field on the sphere
  /// using the analysis and synthesis methods, then computes the analysis and synthesis 
  /// again using the projection field to compute the accuracy of the transform.  Note the
  /// use of the 3-d triangular array to store the spectral coefficients.
  /// \code
  /// #include <iostream>
  /// #include <alloc.h>
  /// #include <sphere.h>
  /// 
  /// using namespace spectral;
  ///
  /// int main()
  /// {	
  ///   int N=32,K=32;
  ///   int threads=2;
  ///   sphere S(N,K,threads);
  ///   real ***g = alloc<real>(K,2*N,4*N);
  ///   real ***g2 = alloc<real>(K,2*N,4*N);
  ///   complex ***sc=alloct<complex>(K,2*N); // triangle
  ///   for(int k=0;k<K;k++)
  ///     for(int i=0;i<2*N;i++)
  ///       for(int j=0;j<4*N;j++)
  ///         g[k][i][j] = (real)(k*j*i)/(real)(K*N*N*8);
  ///   S.analysis(g,sc);
  ///   S.synthesis(sc,g);
  ///   for(int k=0;k<K;k++)
  ///     for(int i=0;i<2*N;i++)		
  ///       for(int j=0;j<4*N;j++)
  ///         g2[k][i][j]=g[k][i][j];	
  ///   S.analysis(g,sc);
  ///   S.synthesis(sc,g);
  ///   real error=0.0;
  ///   for(int k=0;k<K;k++)
  ///     for(int i=0;i<2*N;i++)		
  ///       for(int j=0;j<4*N;j++)
  ///         error=max(error,fabs(g2[k][i][j]-g[k][i][j]));	
  ///   std::cout << " error=" << error << std::endl;
  ///   dealloc<real>(g);
  ///   dealloc<real>(g2);
  ///   dealloc<complex>(sc);
  ///   return 0;
  /// }
  /// \endcode

  class sphere 
  {  
  public:	
    sphere(int N,int K);
    sphere(int N,int K,int threads);
    ~sphere();	
    void analysis(real ***g, complex ***sc);
    void synthesis(complex ***sc, real ***g);
  private:
    /// internal thread object for sphere
    class thread : public fiber
    {
    public:
      thread(int nlat,int ninst,int m0,int ml,int k0,int kl,group *GROUP);
      ~thread();
      void analysis(real ***g,complex ***sc);
      void synthesis(complex ***sc,real ***g);
      void wait();
    private:
      int M,K,N; 
      int M0,ML,K0,KL;
      real *w;
      rfft *RFFT;
      complex ***SC;
      complex ***G;
      real ***g;
      complex ***sc;
      enum method {ANALYSIS,SYNTHESIS} Method; 
      void start();
      group *Group;
      void Analysis();
      void Synthesis();
      void multiply(complex **a,real **b,complex **c,int m,int n,int k);
      real ***P;
      /// recursion relation coefficient array
      struct generator {real c0,c1,c2;} **gen; 
      real ***seed;
      real ****Pmn;
    };
    thread **T;
    int Nodes;
    group *G;
    friend class thread;
  };
}

#endif

// Local Variables:
// mode:C++
// End: