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

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

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

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

20
  size_t ik = 1;
21

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

40
  // 2.0 Compute weights
41

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

  return;
}

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

63
  // 1.0 Initizialization
64

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

68
  size_t iodd = (kn % 2);
69

70
  double zdlx = *pl;
71

72
  // 2.0 Newton iteration
73

74
  for (int jter = 1; jter <= itemax+1; jter++)
75 76 77 78 79 80 81 82 83 84 85 86 87 88
    {
      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
89
void gauaw(size_t kn, double *restrict pl, double *restrict pw)
90 91 92 93 94 95 96
{
  /*
   * 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)
   */
97 98
  double *zfn    = (double *) malloc((kn+1) * (kn+1) * sizeof(double));
  double *zfnlat = (double *) malloc((kn/2+1+1)*sizeof(double));
99 100

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

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

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


120
  // 2.0 Gaussian latitudes and weights
121

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

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

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

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

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

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

148
  // convert to physical latitude
149

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

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

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

  return;
}


167
void gaussianLatitudes(size_t nlats, double *latitudes, double *weights)
168
{
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
  gauaw(nlats, latitudes, weights);
}


bool isGaussianLatitudes(size_t nlats, const double *latitudes)
{
  bool is_gauss_lats = false;

  if (nlats > 2) // check if gaussian
    {
      size_t i;
      double *yv = (double *) malloc(nlats*sizeof(double));
      double *yw = (double *) malloc(nlats*sizeof(double));
      gaussianLatitudes(nlats, yv, yw);
      free(yw);

      for (i = 0; i < nlats; i++)
        yv[i] = asin(yv[i]) / M_PI * 180.0;

      for (i = 0; i < nlats; i++)
        if (fabs(yv[i] - latitudes[i]) > ((yv[0] - yv[1]) / 500.0)) break;

      if (i == nlats) is_gauss_lats = true;

      // check S->N
      if (is_gauss_lats == false)
        {
          for (i = 0; i < nlats; i++)
            if (fabs(yv[i] - latitudes[nlats-i-1]) > ((yv[0] - yv[1]) / 500.0)) break;

          if (i == nlats) is_gauss_lats = true;
        }

      free(yv);
    }

  return is_gauss_lats;
206
}