src/tests/trans/test_transgeneral.cc
//-----------------------------------------------------------------------------
// Routine to compute the spherical harmonics analytically on a grid
// (up to wave number 3)
//
// Author:
// Andreas Mueller *ECMWF*
//
Routine to compute the spherical harmonics analytically on a grid
(up to wave number 3)

Andreas Mueller *ECMWF
Routine to compute the spherical harmonics analytically on a grid ( up to wave number 3 ) Andreas Mueller * ECMWF
void spectral_transform_grid_analytic(const int trc,       
                                      bool trcFT,          
                                      const int n,         
                                      const int m,         
                                      const int imag,      
                                      const Grid grid,     
                                      double rspecg[],     
                                      double rgp[],        
                                      const int ivar_in,   
                                      const int ivar_out)  
{
    int N = (trc + 2) * (trc + 1) / 2;
    for (int jm = 0; jm < 2 * N; jm++) {
        rspecg[jm] = 0.;
    }
    int k = 0;
    for (int jm = 0; jm <= trc; jm++) {
        for (int jn = jm; jn <= trc; jn++) {
            if (jm == m && jn == n) {
                rspecg[2 * k + imag]       = 1.;
                rspecg[2 * k + (1 - imag)] = 0.;
            }
            k++;
        }
    }
    for (int jm = 0; jm < grid.size(); jm++) {
        rgp[jm] = 0.;
    }
    if (StructuredGrid(grid)) {
        StructuredGrid g(grid);
        Grid gridGlobal;
        StructuredGrid gs_global;
        int jlatMin = 0;
        if (trcFT) {
            gridGlobal      = Grid(grid.name());
            gs_global       = StructuredGrid(gridGlobal);
            int nlatsGlobal = gs_global.ny();
            for (int jlat = 0; jlat < nlatsGlobal; jlat++) {
                if (gs_global.y(jlat) > g.y(0)) {
                    jlatMin++;
                }
            }
        }
        int idx = 0;
        for (idx_t j = 0; j < g.ny(); ++j) {
            double lat = g.y(j) * util::Constants::degreesToRadians();
            int ftrc   = trc + 1;
            if (trcFT) {
                ftrc = trans::fourier_truncation(trc, gs_global.nx(jlatMin + j), gs_global.nxmax(), gs_global.ny(), lat,
                                                 RegularGrid(gs_global));
            }
            for (idx_t i = 0; i < g.nx(j); ++i) {
                double lon = g.x(i, j) * util::Constants::degreesToRadians();
                if (ftrc > m) {
                    rgp[idx++] = sphericalharmonics_analytic_point(n, m, imag, lon, lat, ivar_in, ivar_out);
                }
                else {
                    rgp[idx++] = 0.;
                }
            }
        }
    }
    else {
        int idx = 0;
        for (const PointXY& p : grid.xy()) {
            double lon = p.x() * util::Constants::degreesToRadians();
            double lat = p.y() * util::Constants::degreesToRadians();
            rgp[idx++] = sphericalharmonics_analytic_point(n, m, imag, lon, lat, ivar_in, ivar_out);
        }
    }
}
void spectral_transform_grid_analytic ( const int trc , bool trcFT , const int n , const int m , const int imag , const Grid grid , double rspecg [ ] , double rgp [ ] , const int ivar_in , const int ivar_out ) { int N = ( trc + 2 ) * ( trc + 1 ) / 2 ; for ( int jm = 0 ; jm < 2 * N ; jm ++ ) { rspecg [ jm ] = 0. ; } int k = 0 ; for ( int jm = 0 ; jm <= trc ; jm ++ ) { for ( int jn = jm ; jn <= trc ; jn ++ ) { if ( jm == m && jn == n ) { rspecg [ 2 * k + imag ] = 1. ; rspecg [ 2 * k + ( 1 - imag ) ] = 0. ; } k ++ ; } } for ( int jm = 0 ; jm < grid . size ( ) ; jm ++ ) { rgp [ jm ] = 0. ; } if ( StructuredGrid ( grid ) ) { StructuredGrid g ( grid ) ; Grid gridGlobal ; StructuredGrid gs_global ; int jlatMin = 0 ; if ( trcFT ) { gridGlobal = Grid ( grid . name ( ) ) ; gs_global = StructuredGrid ( gridGlobal ) ; int nlatsGlobal = gs_global . ny ( ) ; for ( int jlat = 0 ; jlat < nlatsGlobal ; jlat ++ ) { if ( gs_global . y ( jlat ) > g . y ( 0 ) ) { jlatMin ++ ; } } } int idx = 0 ; for ( idx_t j = 0 ; j < g . ny ( ) ; ++ j ) { double lat = g . y ( j ) * util :: Constants :: degreesToRadians ( ) ; int ftrc = trc + 1 ; if ( trcFT ) { ftrc = trans :: fourier_truncation ( trc , gs_global . nx ( jlatMin + j ) , gs_global . nxmax ( ) , gs_global . ny ( ) , lat , RegularGrid ( gs_global ) ) ; } for ( idx_t i = 0 ; i < g . nx ( j ) ; ++ i ) { double lon = g . x ( i , j ) * util :: Constants :: degreesToRadians ( ) ; if ( ftrc > m ) { rgp [ idx ++ ] = sphericalharmonics_analytic_point ( n , m , imag , lon , lat , ivar_in , ivar_out ) ; } else { rgp [ idx ++ ] = 0. ; } } } } else { int idx = 0 ; for ( const PointXY & p : grid . xy ( ) ) { double lon = p . x ( ) * util :: Constants :: degreesToRadians ( ) ; double lat = p . y ( ) * util :: Constants :: degreesToRadians ( ) ; rgp [ idx ++ ] = sphericalharmonics_analytic_point ( n , m , imag , lon , lat , ivar_in , ivar_out ) ; } } }
Routine to compute the spherical harmonics analytically on a grid
(up to wave number 3)
Routine to compute the spherical harmonics analytically on a grid ( up to wave number 3 )
// truncation (in) // truncation for Fourier transformation (in) // total wave number (implemented so far for n<4) // zonal wave number (implemented so far for m<4, m<n) // 0: test real part, 1: test imaginary part // call with something like Grid("O32") // spectral data, size (trc+1)*trc (out) // resulting grid point data (out) // variable that is set to 1 for wave number n,m. 0: // vorticity, 1: divergence, 2: scalar // variable returned by this function. 0: u, 1: v, 2: scalar /*Log::info() << "j=" << j << " ftrc=" << ftrc << " trc=" << trc << " nx=" << gs_global.nx( jlatMin + j )
                        << " nxmax=" << gs_global.nxmax() << " nlats=" << gs_global.ny() << " lat=" << g.y( j )
                        << " jlatMin=" << jlatMin << std::endl;*/ // compute spherical harmonics: // compute spherical harmonics: