gaussian_latitudes.c 3.63 KB
Newer Older
1
#include <stdio.h>
2
#include <stdlib.h>
3 4
#include <float.h>
#include <math.h>
5

6 7 8
#ifndef  M_SQRT2
#define  M_SQRT2  1.41421356237309504880168872420969808
#endif
Oliver Heidmann's avatar
Oliver Heidmann committed
9

10
static
11
void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
12 13
            double *pw, double *pdxn, double *pxmod)
{
14
  // 1.0 Newton iteration step
15

16 17
  double zdlx = pdx;
  double zdlldn = 0.0;
18

19
  size_t ik = 1;
20

21
  if (kflag == 0)
22
    {
23 24
      double zdlk = 0.5*pfn[0];
      for (size_t jn = 2-kodd; jn <= kn; jn += 2)
25
	{
26
	  // normalised ordinary Legendre polynomial == \overbar{p_n}^0
27
	  zdlk   = zdlk + pfn[ik]*cos((double)(jn)*zdlx);
28
	  // normalised derivative == d/d\theta(\overbar{p_n}^0)
29 30 31
	  zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
	  ik++;
	}
32
      // Newton method
33
      double zdlmod = -(zdlk/zdlldn);
34
      double zdlxn = zdlx + zdlmod;
35 36 37 38
      *pdxn = zdlxn;
      *pxmod = zdlmod;
    }

39
  // 2.0 Compute weights
40

41
  if ( kflag == 1 )
42
    {
43
      for (size_t jn = 2-kodd; jn <= kn; jn += 2)
44
	{
45
	  // normalised derivative
46 47 48 49 50 51 52 53 54 55
	  zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
	  ik++;
	}
      *pw = (double)(2*kn+1)/(zdlldn*zdlldn);
    }

  return;
}

static
56
void gawl(double *pfn, double *pl, double *pw, size_t kn)
57
{
58 59 60
  double pmod = 0.0;
  double zw = 0.0;
  double zdlxn = 0.0;
61

62
  // 1.0 Initizialization
63

64 65
  int iflag  =  0;
  int itemax = 20;
66

67
  size_t iodd = (kn % 2);
68

69
  double zdlx = *pl;
70

71
  // 2.0 Newton iteration
72

73
  for (int jter = 1; jter <= itemax+1; jter++)
74 75 76 77 78 79 80 81 82 83 84 85 86 87
    {
      cpledn(kn, iodd, pfn, zdlx, iflag, &zw, &zdlxn, &pmod);
      zdlx = zdlxn;
      if (iflag == 1) break;
      if (fabs(pmod) <= DBL_EPSILON*1000.0) iflag = 1;
    }

  *pl = zdlxn;
  *pw = zw;

  return;
}

static
88
void gauaw(size_t kn, double *restrict pl, double *restrict pw)
89 90 91 92 93 94 95
{
  /*
   * 1.0 Initialize Fourier coefficients for ordinary Legendre polynomials
   *
   * Belousov, Swarztrauber, and ECHAM use zfn(0,0) = sqrt(2)
   * IFS normalisation chosen to be 0.5*Integral(Pnm**2) = 1 (zfn(0,0) = 2.0)
   */
96 97
  double *zfn    = (double *) malloc((kn+1) * (kn+1) * sizeof(double));
  double *zfnlat = (double *) malloc((kn/2+1+1)*sizeof(double));
98 99

  zfn[0] = M_SQRT2;
100
  for (size_t jn = 1; jn <= kn; jn++)
101
    {
102
      double zfnn = zfn[0];
103
      for (size_t jgl = 1; jgl <= jn; jgl++)
104
	{
105
	  zfnn *= sqrt(1.0-0.25/((double)(jgl*jgl)));
106 107 108 109
	}

      zfn[jn*(kn+1)+jn] = zfnn;

110
      size_t iodd = jn % 2;
111
      for (size_t jgl = 2; jgl <= jn-iodd; jgl += 2)
112 113
	{
	  zfn[jn*(kn+1)+jn-jgl] = zfn[jn*(kn+1)+jn-jgl+2]
114
	    * ((double)((jgl-1)*(2*jn-jgl+2))) / ((double)(jgl*(2*jn-jgl+1)));
115 116 117 118
	}
    }


119
  // 2.0 Gaussian latitudes and weights
120

121 122
  size_t iodd = kn % 2;
  size_t ik = iodd;
123
  for (size_t jgl = iodd; jgl <= kn; jgl += 2)
124 125 126
    {
      zfnlat[ik] = zfn[kn*(kn+1)+jgl];
      ik++;
127
    }
128

129
  // 2.1 Find first approximation of the roots of the Legendre polynomial of degree kn
130

131
  size_t ins2 = kn/2+(kn % 2);
132

133
  for (size_t jgl = 1; jgl <= ins2; jgl++)
134
    {
135 136
      double z = ((double)(4*jgl-1)) * M_PI / ((double)(4*kn+2));
      pl[jgl-1] = z + 1.0 / (tan(z)*((double)(8*kn*kn)));
137 138
    }

139
  // 2.2 Computes roots and weights for transformed theta
140

141
  for (size_t jgl = ins2; jgl >= 1 ; jgl--)
142
    {
143
      size_t jglm1 = jgl-1;
Thomas Jahns's avatar
Thomas Jahns committed
144
      gawl(zfnlat, &(pl[jglm1]), &(pw[jglm1]), kn);
145 146
    }

147
  // convert to physical latitude
148

149
  for (size_t jgl = 0; jgl < ins2; jgl++) pl[jgl] = cos(pl[jgl]);
150

151
  for (size_t jgl = 1; jgl <= kn/2; jgl++)
152
    {
153
      size_t jglm1 = jgl-1;
154 155
      size_t isym  = kn-jgl;
      pl[isym] = -pl[jglm1];
156 157 158
      pw[isym] =  pw[jglm1];
    }

159 160
  free(zfnlat);
  free(zfn);
161 162 163 164 165

  return;
}


166
void gaussianLatitudes(double *latitudes, double *weights, size_t nlat)
167
{
168
  gauaw(nlat, latitudes, weights);
169
}