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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
#include "dmemory.h"

8
9
10
#ifndef  M_PI
#define  M_PI        3.14159265358979323846  /* pi */
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
11
12
13
14
15
16

#ifndef  M_SQRT2
#define  M_SQRT2     1.41421356237309504880
#endif


17
static
18
void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag, 
19
20
21
22
23
24
25
26
            double *pw, double *pdxn, double *pxmod)
{
  double zdlk;
  double zdlldn;
  double zdlx;
  double zdlmod;
  double zdlxn;

27
  size_t ik, jn;
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

  /* 1.0 Newton iteration step */

  zdlx = pdx;
  zdlk = 0.0;
  if (kodd == 0) 
    {
      zdlk = 0.5*pfn[0];
    }
  zdlxn  = 0.0;
  zdlldn = 0.0;

  ik = 1;

  if (kflag == 0) 
    {
44
      for(size_t jn = 2-kodd; jn <= kn; jn += 2) 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
	{
	  /* normalised ordinary Legendre polynomial == \overbar{p_n}^0 */
	  zdlk   = zdlk + pfn[ik]*cos((double)(jn)*zdlx);
	  /* normalised derivative == d/d\theta(\overbar{p_n}^0) */
	  zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
	  ik++;
	}
      /* Newton method */
      zdlmod = -(zdlk/zdlldn);
      zdlxn = zdlx + zdlmod;
      *pdxn = zdlxn;
      *pxmod = zdlmod;
    }

  /* 2.0 Compute weights */

  if (kflag == 1) 
    {
      for(jn = 2-kodd; jn <= kn; jn += 2) 
	{
	  /* normalised derivative */
	  zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
	  ik++;
	}
      *pw = (double)(2*kn+1)/(zdlldn*zdlldn);
    }

  return;
}

static
76
void gawl(double *pfn, double *pl, double *pw, size_t kn)
77
78
79
80
81
82
83
84
85
86
87
88
89
{
  double pmod = 0;
  int iflag;
  int itemax;
  double zw = 0;
  double zdlx;
  double zdlxn = 0;

  /* 1.0 Initizialization */

  iflag  =  0;
  itemax = 20;

90
  size_t iodd   = (kn % 2);
91
92
93
94
95

  zdlx   =  *pl;

  /* 2.0 Newton iteration */

Thomas Jahns's avatar
Thomas Jahns committed
96
  for (int jter = 1; jter <= itemax+1; jter++)
97
98
99
100
101
102
103
104
105
106
107
108
109
110
    {
      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
111
void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
112
113
114
115
116
117
118
119
120
121
122
{
  /*
   * 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)
   */
  double *zfn, *zfnlat;

  double z, zfnn;

123
124
  zfn    = (double *) malloc((kn+1) * (kn+1) * sizeof(double));
  zfnlat = (double *) malloc((kn/2+1+1)*sizeof(double));
125
126

  zfn[0] = M_SQRT2;
127
  for (size_t jn = 1; jn <= kn; jn++)
128
129
    {
      zfnn = zfn[0];
130
      for (size_t jgl = 1; jgl <= jn; jgl++)
131
132
133
134
135
136
	{
	  zfnn *= sqrt(1.0-0.25/((double)(jgl*jgl))); 
	}

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

137
138
      size_t iodd = jn % 2;
      for (size_t jgl = 2; jgl <= jn-iodd; jgl += 2) 
139
140
141
142
143
144
145
146
147
	{
	  zfn[jn*(kn+1)+jn-jgl] = zfn[jn*(kn+1)+jn-jgl+2]
	    *((double)((jgl-1)*(2*jn-jgl+2)))/((double)(jgl*(2*jn-jgl+1)));
	}
    }


  /* 2.0 Gaussian latitudes and weights */

148
149
150
  size_t iodd = kn % 2;
  size_t ik = iodd;
  for (size_t jgl = iodd; jgl <= kn; jgl += 2)
151
152
153
154
155
156
157
158
159
160
    {
      zfnlat[ik] = zfn[kn*(kn+1)+jgl];
      ik++;
    } 

  /*
   * 2.1 Find first approximation of the roots of the
   *     Legendre polynomial of degree kn.
   */

161
  size_t ins2 = kn/2+(kn % 2);
162

163
  for (size_t jgl = 1; jgl <= ins2; jgl++) 
164
165
166
167
168
169
170
    {
      z = ((double)(4*jgl-1))*M_PI/((double)(4*kn+2)); 
      pl[jgl-1] = z+1.0/(tan(z)*((double)(8*kn*kn)));
    }

  /* 2.2 Computes roots and weights for transformed theta */

171
  for (size_t jgl = ins2; jgl >= 1 ; jgl--) 
172
    {
173
      size_t jglm1 = jgl-1;
Thomas Jahns's avatar
Thomas Jahns committed
174
      gawl(zfnlat, &(pl[jglm1]), &(pw[jglm1]), kn);
175
176
177
178
    }

  /* convert to physical latitude */

179
  for (size_t jgl = 0; jgl < ins2; jgl++) 
180
181
182
183
    {
      pl[jgl] = cos(pl[jgl]);
    }

184
  for (size_t jgl = 1; jgl <= kn/2; jgl++) 
185
    {
186
187
      size_t jglm1 = jgl-1;
      size_t isym =  kn-jgl;
188
189
190
191
192
193
194
195
196
197
      pl[isym] =  -pl[jglm1];
      pw[isym] =  pw[jglm1];
    }

  free(zfnlat);
  free(zfn);

  return;
}

198
#if 0
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
static
void gauaw_old(double *pa, double *pw, int nlat)
{
  /*
   * Compute Gaussian latitudes.  On return pa contains the
   * sine of the latitudes starting closest to the north pole and going
   * toward the south
   *
   */

  const int itemax = 20;

  int isym, iter, ins2, jn, j;
  double za, zw, zan;
  double z, zk, zkm1, zkm2, zx, zxn, zldn, zmod;

  /*
   * Perform the Newton loop
   * Find 0 of Legendre polynomial with Newton loop
   */

  ins2 = nlat/2 + nlat%2;

  for ( j = 0; j < ins2; j++ )
    {
      z = (double) (4*(j+1)-1)*M_PI / (double) (4*nlat+2);
      pa[j] = cos(z + 1.0/(tan(z)*(double)(8*nlat*nlat)));
    }

  for ( j = 0; j < ins2; j++ )
    {

      za = pa[j];

      iter = 0;
      do
	{
	  iter++;
	  zk = 0.0;

	  /* Newton iteration step */

	  zkm2 = 1.0;
	  zkm1 = za;
	  zx = za;
	  for ( jn = 2; jn <= nlat; jn++ )
	    {
	      zk = ((double) (2*jn-1)*zx*zkm1-(double)(jn-1)*zkm2) / (double)(jn);
	      zkm2 = zkm1;
	      zkm1 = zk;
	    }
	  zkm1 = zkm2;
	  zldn = ((double) (nlat)*(zkm1-zx*zk)) / (1.-zx*zx);
	  zmod = -zk/zldn;
	  zxn = zx+zmod;
	  zan = zxn;

	  /* computes weight */

	  zkm2 = 1.0;
	  zkm1 = zxn;
	  zx = zxn;
	  for ( jn = 2; jn <= nlat; jn++ )
	    {
	      zk = ((double) (2*jn-1)*zx*zkm1-(double)(jn-1)*zkm2) / (double) (jn);
	      zkm2 = zkm1;
	      zkm1 = zk;
	    }
	  zkm1 = zkm2;
	  zw = (1.0-zx*zx) / ((double) (nlat*nlat)*zkm1*zkm1);
	  za = zan;
	}
      while ( iter <= itemax && fabs(zmod) >= DBL_EPSILON );

      pa[j] = zan;
      pw[j] = 2.0*zw;
    }

#if defined (SX)
#pragma vdir nodep
#endif
  for (j = 0; j < nlat/2; j++)
    {
      isym = nlat-(j+1);
      pa[isym] = -pa[j];
      pw[isym] =  pw[j];
    }

  return;
}
289
#endif
290

291
void gaussaw(double *restrict pa, double *restrict pw, size_t nlat)
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
{
  //gauaw_old(pa, pw, nlat);
  gauaw(nlat, pa, pw);
}

/*
#define NGL  48

int main (int rgc, char *argv[])
{
  int ngl = NGL;
  double plo[NGL], pwo[NGL];
  double pl[NGL], pw[NGL];

  int i;

  gauaw(ngl, pl, pw);
  gauaw_old(plo, pwo, ngl);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310
311
312
313
314
  for (i = 0; i < ngl; i++)
    {
      pl[i]  = asin(pl[i])/M_PI*180.0;
      plo[i] = asin(plo[i])/M_PI*180.0;
    }
315
316
317

  for (i = 0; i < ngl; i++)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318
      fprintf(stderr, "%4d%25.18f%25.18f%25.18f%25.18f\n", i+1, pl[i], pw[i], pl[i]-plo[i], pw[i]-pwo[i]);
319
320
321
322
323
    }

  return 0;
}
*/
324
325
326
327
328
329
330
331
332
/*
 * Local Variables:
 * c-file-style: "Java"
 * c-basic-offset: 2
 * indent-tabs-mode: nil
 * show-trailing-whitespace: t
 * require-trailing-newline: t
 * End:
 */