ucar.unidata.util.GaussianLatitudes Maven / Gradle / Ivy
The newest version!
/*
* Copyright (c) 1998-2018 University Corporation for Atmospheric Research/Unidata
* See LICENSE for license information.
*/
package ucar.unidata.util;
import ucar.nc2.util.HashMapLRU;
import javax.annotation.concurrent.Immutable;
/**
* Calculate Gaussian Latitudes by finding the roots of the ordinary Legendre polynomial of degree
* NLAT using Newton's iteration method.
* Moderately expensive to compute, so keep LRU cache.
*
* @author caron port to Java
* @author Amy Solomon, 28 Jan 1991, via http://dss.ucar.edu/libraries/gridinterps/gaus-lats.f
*/
@Immutable
public class GaussianLatitudes {
private static final double XLIM = 1.0E-7;
private static final HashMapLRU hash = new HashMapLRU<>(10, 20);
public static GaussianLatitudes factory(int nlats) {
GaussianLatitudes glats = hash.get(nlats);
if (glats == null) {
glats = new GaussianLatitudes(nlats);
hash.put(nlats, glats);
}
return glats;
}
public final double[] gaussw; // the Gaussian weights
public final double[] latd; // the latitudes in degrees
/**
* Constructor
*
* @param nlat the total number of latitudes from pole to pole (degree of the polynomial)
*/
private GaussianLatitudes(int nlat) {
if (nlat == 0)
throw new IllegalArgumentException("nlats may not be zero");
// the number of latitudes between pole and equator
int nzero = nlat / 2;
// all these have size nlat
// cos(colatitude) or sin(latitude)
double[] cosc = new double[nlat];
// the colatitudes in radians
double[] colat = new double[nlat];
gaussw = new double[nlat];
latd = new double[nlat];
// sin(colatitude) or cos(latitude)
double[] sinc = new double[nlat];
// Gaussian weight over sin**2(colatitude)
double[] wos2 = new double[nlat];
/*
* set first guess for cos(colat)
* PI = 3.141592653589793
* DO 10 I=1,NZERO
* COSC(I) = SIN( (I-0.5)*PI/NLAT + PI*0.5 )
* 10 CONTINUE
*/
for (int i = 0; i < nzero; i++)
cosc[i] = Math.sin((i + 0.5) * Math.PI / nlat + Math.PI / 2);
/*
* constants for determining the derivative of the polynomial
* FI = NLAT
* FI1 = FI+1.0
* A = FI*FI1 / SQRT(4.0*FI1*FI1-1.0)
* B = FI1*FI / SQRT(4.0*FI*FI-1.0)
*/
double FI = (double) nlat;
double FI1 = nlat + 1.0;
double A = FI * FI1 / Math.sqrt(4.0 * FI1 * FI1 - 1.0);
double B = FI * FI1 / Math.sqrt(4.0 * FI * FI - 1.0);
// loop over latitudes, iterating the search for each root
for (int i = 0; i < nzero; i++) {
int countIterations = 0;
while (true) {
// -determine the value of the ordinary Legendre polynomial for the current guess root
/* 30 CALL LGORD( G, COSC(I), NLAT ) */
double G = lgord(cosc[i], nlat);
/*
* -determine the derivative of the polynomial at this point
* CALL LGORD( GM, COSC(I), NLAT-1 )
* CALL LGORD( GP, COSC(I), NLAT+1 )
* GT = (COSC(I)*COSC(I)-1.0) / (A*GP-B*GM)
*/
double GM = lgord(cosc[i], nlat - 1);
double GP = lgord(cosc[i], nlat + 1);
double denom = A * GP - B * GM;
double GT = (denom == 0.0) ? 0.0 : (cosc[i] * cosc[i] - 1.0) / denom;
/*
* -update the estimate of the root
* DELTA = G*GT
* COSC(I) = COSC(I) - DELTA
*/
double delta = G * GT;
cosc[i] -= delta;
/*
* -if convergence criterion has not been met, keep trying
* J = J+1
* IF( ABS(DELTA).GT.XLIM ) GO TO 30
*/
countIterations++;
if (Math.abs(delta) <= XLIM)
break;
}
/*
* determine the Gaussian weights
* C = 2.0 *( 1.0-COSC(I)*COSC(I) )
* CALL LGORD( D, COSC(I), NLAT-1 )
* D = D*D*FI*FI
* GWT(I) = C *( FI-0.5 ) / D
*/
double C = 2.0 * (1.0 - cosc[i] * cosc[i]);
double D = lgord(cosc[i], nlat - 1);
D = D * D * FI * FI;
gaussw[i] = (D == 0.0) ? 0.0 : C * (FI - 0.5) / D;
}
/*
* C -determine the colatitudes and sin(colat) and weights over sin**2
* DO 50 I=1,NZERO
* COLAT(I)= ACOS(COSC(I))
* SINC(I) = SIN(COLAT(I))
* WOS2(I) = GWT(I) /( SINC(I)*SINC(I) )
* 50 CONTINUE
*/
for (int i = 0; i < nzero; i++) {
colat[i] = Math.acos(cosc[i]);
sinc[i] = Math.sin(colat[i]);
wos2[i] = gaussw[i] / (sinc[i] * sinc[i]);
}
/*
* if NLAT is odd, set values at the equator
* IF( MOD(NLAT,2) .NE. 0 ) THEN
* I = NZERO+1
* COSC(I) = 0.0
* C = 2.0
* CALL LGORD( D, COSC(I), NLAT-1 )
* D = D*D*FI*FI
* GWT(I) = C *( FI-0.5 ) / D
* COLAT(I)= PI*0.5
* SINC(I) = 1.0
* WOS2(I) = GWT(I)
* END IF
*/
int next = nzero;
if (nlat % 2 != 0) {
cosc[next] = 0.0;
double C = 2.0;
double D = lgord(cosc[next], nlat - 1);
D = D * D * FI * FI;
gaussw[next] = (D == 0.0) ? 0.0 : C * (FI - 0.5) / D;
colat[next] = Math.PI * 0.5;
sinc[next] = 1.0;
wos2[next] = gaussw[next];
next++;
}
/*
* determine the southern hemisphere values by symmetry
* DO 60 I=NLAT-NZERO+1,NLAT
* COSC(I) =-COSC(NLAT+1-I)
* GWT(I) = GWT(NLAT+1-I)
* COLAT(I)= PI-COLAT(NLAT+1-I)
* SINC(I) = SINC(NLAT+1-I)
* WOS2(I) = WOS2(NLAT+1-I)
* 60 CONTINUE
*/
for (int i = next; i < nlat; i++) {
cosc[i] = -cosc[nlat - i - 1];
gaussw[i] = gaussw[nlat - i - 1];
colat[i] = Math.PI - colat[nlat - i - 1];
sinc[i] = sinc[nlat - i - 1];
wos2[i] = wos2[nlat - i - 1];
}
// the lats in degrees
for (int i = 0; i < nlat; i++)
latd[i] = Math.toDegrees(Math.PI / 2 - colat[i]);
}
private double lgord(double cosc, int nlat) {
/*
* COLAT = ACOS(COSC)
* C1 = SQRT(2.0)
* DO 20 K=1,N
* C1 = C1 * SQRT( 1.0 - 1.0/(4*K*K) )
* 20 CONTINUE
*/
double colat = Math.acos(cosc);
double c = Math.sqrt(2.0);
for (int k = 1; k <= nlat; k++)
c *= Math.sqrt(1.0 - 1.0 / (4 * k * k));
/*
* FN = N
* ANG= FN * COLAT
* S1 = 0.0
* C4 = 1.0
* A =-1.0
* B = 0.0
* DO 30 K=0,N,2
* IF (K.EQ.N) C4 = 0.5 * C4
* S1 = S1 + C4 * COS(ANG)
* A = A + 2.0
* B = B + 1.0
* FK = K
* ANG= COLAT * (FN-FK-2.0)
* C4 = ( A * (FN-B+1.0) / ( B * (FN+FN-A) ) ) * C4
* 30 CONTINUE
* F = S1 * C1
*/
double FN = (double) nlat;
double ANG = FN * colat;
double S1 = 0.0;
double C4 = 1.0;
double A = -1.0;
double B = 0.0;
for (int k = 0; k <= nlat; k += 2) {
if (k == nlat)
C4 = 0.5 * C4;
S1 = S1 + C4 * Math.cos(ANG);
A = A + 2.0;
B = B + 1.0;
double FK = (double) k;
ANG = colat * (FN - FK - 2.0);
C4 = (A * (FN - B + 1.0) / (B * (FN + FN - A))) * C4;
}
return S1 * c;
}
}
/*
*
* C ********** WARNING *************
* C This routine may not converge using 32-bit arithmetic
* C
* C routines from Amy Solomon, 28 Jan 1991.
* C
* C LGGAUS finds the Gaussian latitudes by finding the roots of the
* C ordinary Legendre polynomial of degree NLAT using Newton's
* C iteration method.
* C
* C On entry:
* C NLAT - the number of latitudes (degree of the polynomial)
* C
* C On exit: for each Gaussian latitude
* C COSC - cos(colatitude) or sin(latitude)
* C GWT - the Gaussian weights
* C SINC - sin(colatitude) or cos(latitude)
* C COLAT - the colatitudes in radians
* C WOS2 - Gaussian weight over sin**2(colatitude)
* C
* DIMENSION COSC(180), GWT(180), SINC(180), COLAT(180)
* + , WOS2(180)
* C
* C-----------------------------------------------------------------------
* C
* C -convergence criterion for iteration of cos latitude
* XLIM = 1.0E-7
* C
* C -the number of zeros between pole and equator
* NZERO = NLAT/2
* C
* C -set first guess for cos(colat)
* PI = 3.141592653589793
* DO 10 I=1,NZERO
* COSC(I) = SIN( (I-0.5)*PI/NLAT + PI*0.5 )
* 10 CONTINUE
* C
* C -constants for determining the derivative of the polynomial
* FI = NLAT
* FI1 = FI+1.0
* A = FI*FI1 / SQRT(4.0*FI1*FI1-1.0)
* B = FI1*FI / SQRT(4.0*FI*FI-1.0)
* C
* C -loop over latitudes, iterating the search for each root
* DO 40 I=1,NZERO
* J=0
* C
* C -determine the value of the ordinary Legendre polynomial for
* C -the current guess root
* 30 CALL LGORD( G, COSC(I), NLAT )
* C
* C -determine the derivative of the polynomial at this point
* CALL LGORD( GM, COSC(I), NLAT-1 )
* CALL LGORD( GP, COSC(I), NLAT+1 )
* GT = (COSC(I)*COSC(I)-1.0) / (A*GP-B*GM)
* C
* C -update the estimate of the root
* DELTA = G*GT
* COSC(I) = COSC(I) - DELTA
* C
* C -if convergence criterion has not been met, keep trying
* J = J+1
* IF( ABS(DELTA).GT.XLIM ) GO TO 30
* C PRINT*,' LAT NO.',I,J,' ITERATIONS'
* C
* C -determine the Gaussian weights
* C = 2.0 *( 1.0-COSC(I)*COSC(I) )
* CALL LGORD( D, COSC(I), NLAT-1 )
* D = D*D*FI*FI
* GWT(I) = C *( FI-0.5 ) / D
* 40 CONTINUE
* C
* C -determine the colatitudes and sin(colat) and weights over sin**2
* DO 50 I=1,NZERO
* COLAT(I)= ACOS(COSC(I))
* SINC(I) = SIN(COLAT(I))
* WOS2(I) = GWT(I) /( SINC(I)*SINC(I) )
* 50 CONTINUE
* C
* C -if NLAT is odd, set values at the equator
* IF( MOD(NLAT,2) .NE. 0 ) THEN
* I = NZERO+1
* COSC(I) = 0.0
* C = 2.0
* CALL LGORD( D, COSC(I), NLAT-1 )
* D = D*D*FI*FI
* GWT(I) = C *( FI-0.5 ) / D
* COLAT(I)= PI*0.5
* SINC(I) = 1.0
* WOS2(I) = GWT(I)
* END IF
* C
* C -determine the southern hemisphere values by symmetry
* DO 60 I=NLAT-NZERO+1,NLAT
* COSC(I) =-COSC(NLAT+1-I)
* GWT(I) = GWT(NLAT+1-I)
* COLAT(I)= PI-COLAT(NLAT+1-I)
* SINC(I) = SINC(NLAT+1-I)
* WOS2(I) = WOS2(NLAT+1-I)
* 60 CONTINUE
* c PRINT*,'NLAT=',NLAT
* c PRINT*,'COLATS'
* c PRINT 101,(I,COLAT(I),COLAT(I)*180./PI,I=1,NLAT)
* 101 FORMAT(1X,I3,F16.12,2X,F8.2)
* 102 FORMAT(1X,I3,F16.12,2X,F16.12)
* c PRINT*,'COS(COLAT), SIN(COLAT)'
* c PRINT 102,(I,COSC(I),SINC(I),I=1,NLAT)
* c PRINT*,'WEIGHT, GWT/COS**2'
* c PRINT 102,(I,GWT(I),WOS2(I),I=1,NLAT)
* C
* RETURN
* END
* SUBROUTINE LGORD( F, COSC, N )
* C
* C LGORD calculates the value of an ordinary Legendre polynomial at a
* C latitude.
* C
* C On entry:
* C COSC - cos(colatitude)
* C N - the degree of the polynomial
* C
* C On exit:
* C F - the value of the Legendre polynomial of degree N at
* C latitude asin(COSC)
* C
* C------------------------------------------------------------------------
* C
* C -determine the colatitude
* COLAT = ACOS(COSC)
* C
* C1 = SQRT(2.0)
* DO 20 K=1,N
* C1 = C1 * SQRT( 1.0 - 1.0/(4*K*K) )
* 20 CONTINUE
* C
* FN = N
* ANG= FN * COLAT
* S1 = 0.0
* C4 = 1.0
* A =-1.0
* B = 0.0
* DO 30 K=0,N,2
* IF (K.EQ.N) C4 = 0.5 * C4
* S1 = S1 + C4 * COS(ANG)
* A = A + 2.0
* B = B + 1.0
* FK = K
* ANG= COLAT * (FN-FK-2.0)
* C4 = ( A * (FN-B+1.0) / ( B * (FN+FN-A) ) ) * C4
* 30 CONTINUE
* F = S1 * C1
* C
* RETURN
* END
*/
© 2015 - 2025 Weber Informatics LLC | Privacy Policy