qpms/tests/gsltest.c

   1 #include <stdio.h>
   2 #include <stddef.h>
   3 #include <stdlib.h>
   4 #include <gsl/gsl_sf_legendre.h>
   5
   6
   7 int main(int argc, char **argv) {
   8         int lmax;
   9         int scanned;
  10         double x;
  11         while (EOF != (scanned = scanf("%d %lf", &lmax, &x))) {
  12                 if (scanned != 2) continue;
  13                 size_t as = gsl_sf_legendre_array_n(lmax);
  14                 double *r_none = calloc(as, sizeof(double));
  15                 double *d_none = calloc(as, sizeof(double));
  16                 double *r_schmidt = calloc(as, sizeof(double));
  17                 double *d_schmidt = calloc(as, sizeof(double));
  18                 double *r_spharm = calloc(as, sizeof(double));
  19                 double *d_spharm = calloc(as, sizeof(double));
  20                 double *r_full = calloc(as, sizeof(double));
  21                 double *d_full = calloc(as, sizeof(double));
  22                 if( 0
  23                                 || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_NONE,lmax,x,-1,r_none,d_none)
  24                                 || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_SCHMIDT,lmax,x,-1,r_schmidt,d_schmidt)
  25                                 || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_SPHARM,lmax,x,-1,r_spharm,d_spharm)
  26                                 || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_FULL,lmax,x,-1,r_full,d_full)
  27                   ) fprintf(stderr, "Something gone wrong for lmax = %d, x = %.15e!\n", lmax, x);
  28                 for (int l = 0; l <= lmax; ++l)
  29                         for (int m = 0; m <= l; ++m){
  30                                 size_t i = gsl_sf_legendre_array_index(l,m);
  31                                 printf("P(%d,%d)\t%.16e\t%.16e\t%.16e\t%.16e\n", l, m,
  32                                                 r_none[i], r_schmidt[i], r_spharm[i], r_full[i]);
  33                                 printf("dP(%d,%d)\t%.16e\t%.16e\t%.16e\t%.16e\n", l, m,
  34                                                 d_none[i], d_schmidt[i], d_spharm[i], d_full[i]);
  35
  36                         }
  37
  38
  39
  40
  41                 free(r_none);
  42                 free(d_none);
  43                 free(r_schmidt);
  44                 free(d_schmidt);
  45                 free(r_spharm);
  46                 free(d_spharm);
  47                 free(r_full);
  48                 free(d_full);
  49         }
  50         return 0;
  51 }
  52
  53
  54