//// 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: