From: curt Date: Tue, 28 Mar 2000 15:20:20 +0000 (+0000) Subject: Optimizations by Norman Vine: X-Git-Url: https://git.mxchange.org/?a=commitdiff_plain;h=3dde4113e7d0b0bc970d56bf42e97f2c1c8f0f98;p=simgear.git Optimizations by Norman Vine: Classic space vs time seemed worth it in that we get a ~3 fold speedup for ~5% space increase here. Also pow() is an expensive Fortran to C translation in just about all the old government code I see :)) --- diff --git a/simgear/magvar/magvar.cxx b/simgear/magvar/magvar.cxx index 18ae58c1..319dca13 100644 --- a/simgear/magvar/magvar.cxx +++ b/simgear/magvar/magvar.cxx @@ -16,6 +16,11 @@ // Released under GPL 3/26/00 EAW // Adaptions and modifications for the SimGear project 3/27/2000 CLO // +// Removed all pow() calls and made static roots[][] arrays to +// save many sqrt() calls on subsequent invocations +// left old code as SGMagVarOrig() for testing purposes +// 3/28/2000 Norman Vine -- nhv@yahoo.com +// // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License as // published by the Free Software Foundation; either version 2 of the @@ -126,6 +131,8 @@ static double hnm[13][13]; static double sm[13]; static double cm[13]; +static double root[13]; +static double roots[13][13][2]; /* Convert date to Julian day 1950-2049 */ unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd ) @@ -156,9 +163,10 @@ double rad_to_deg( double rad ) } -/* return variation (in degrees) given geodetic latitude (radians), longitude -(radians) ,height (km) and (Julian) date -N and E lat and long are positive, S and W negative +/* + * return variation (in radians) given geodetic latitude (radians), + * longitude(radians), height (km) and (Julian) date + * N and E lat and long are positive, S and W negative */ double SGMagVar( double lat, double lon, double h, long dat, double* field ) @@ -168,6 +176,156 @@ double SGMagVar( double lat, double lon, double h, long dat, double* field ) /* reference dates */ long date0_wmm2000 = yymmdd_to_julian_days(0,1,1); + double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z; + double sinpsi, cospsi, inv_s; + + static int been_here = 0; + + double sinlat = sin(lat); + double coslat = cos(lat); + + /* convert to geocentric coords: */ + // sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0)); + sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat); + /* sr is effective radius */ + theta = atan2(coslat * (h*sr + a*a), + sinlat * (h*sr + b*b)); + /* theta is geocentric co-latitude */ + + r = h*h + 2.0*h * sr + + (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) / + (a*a - (a*a - b*b) * sinlat*sinlat ); + + r = sqrt(r); + + /* r is geocentric radial distance */ + c = cos(theta); + s = sin(theta); + inv_s = 1.0 / s; + + /*zero out arrays */ + for ( n = 0; n <= nmax; n++ ) { + for ( m = 0; m <= n; m++ ) { + P[n][m] = 0; + DP[n][m] = 0; + } + } + + /* diagonal elements */ + P[0][0] = 1; + P[1][1] = s; + DP[0][0] = 0; + DP[1][1] = c; + P[1][0] = c ; + DP[1][0] = -s; + + // these values will not change for subsequent function calls + if( !been_here ) { + for ( n = 2; n <= nmax; n++ ) { + root[n] = sqrt((2.0*n-1) / (2.0*n)); + } + + for ( m = 0; m <= nmax; m++ ) { + double mm = m*m; + for ( n = max(m + 1, 2); n <= nmax; n++ ) { + roots[m][n][0] = sqrt((n-1)*(n-1) - mm); + roots[m][n][1] = 1.0 / sqrt( n*n - mm); + } + } + been_here = 1; + } + + for ( n=2; n <= nmax; n++ ) { + // double root = sqrt((2.0*n-1) / (2.0*n)); + P[n][n] = P[n-1][n-1] * s * root[n]; + DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) * + root[n]; + } + + /* lower triangle */ + for ( m = 0; m <= nmax; m++ ) { + // double mm = m*m; + for ( n = max(m + 1, 2); n <= nmax; n++ ) { + // double root1 = sqrt((n-1)*(n-1) - mm); + // double root2 = 1.0 / sqrt( n*n - mm); + P[n][m] = (P[n-1][m] * c * (2.0*n-1) - + P[n-2][m] * roots[m][n][0]) * + roots[m][n][1]; + + DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) * + (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) * + roots[m][n][1]; + } + } + + /* compute gnm, hnm at dat */ + /* WMM2000 */ + yearfrac = (dat - date0_wmm2000) / 365.25; + for ( n = 1; n <= nmax; n++ ) { + for ( m = 0; m <= nmax; m++ ) { + gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m]; + hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m]; + } + } + + /* compute sm (sin(m lon) and cm (cos(m lon)) */ + for ( m = 0; m <= nmax; m++ ) { + sm[m] = sin(m * lon); + cm[m] = cos(m * lon); + } + + /* compute B fields */ + B_r = 0.0; + B_theta = 0.0; + B_phi = 0.0; + fn_0 = r_0/r; + fn = fn_0 * fn_0; + + for ( n = 1; n <= nmax; n++ ) { + double c1_n=0; + double c2_n=0; + double c3_n=0; + for ( m = 0; m <= n; m++ ) { + double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]); + c1_n=c1_n + tmp * P[n][m]; + c2_n=c2_n + tmp * DP[n][m]; + c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m]; + } + // fn=pow(r_0/r,n+2.0); + fn *= fn_0; + B_r = B_r + (n + 1) * c1_n * fn; + B_theta = B_theta - c2_n * fn; + B_phi = B_phi + c3_n * fn * inv_s; + } + + /* Find geodetic field components: */ + psi = theta - ((pi / 2.0) - lat); + sinpsi = sin(psi); + cospsi = cos(psi); + X = -B_theta * cospsi - B_r * sinpsi; + Y = B_phi; + Z = B_theta * sinpsi - B_r * cospsi; + + field[0]=B_r; + field[1]=B_theta; + field[2]=B_phi; + field[3]=X; + field[4]=Y; + field[5]=Z; /* output fields */ + + /* find variation, leave in radians! */ + return atan2(Y, X); /* E is positive */ +} + + +#ifdef TEST_NHV_HACKS +double SGMagVarOrig( double lat, double lon, double h, long dat, double* field ) +{ + /* output field B_r,B_th,B_phi,B_x,B_y,B_z */ + int n,m; + /* reference dates */ + long date0_wmm2000 = yymmdd_to_julian_days(0,1,1); + double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z; /* convert to geocentric coords: */ @@ -183,11 +341,11 @@ double SGMagVar( double lat, double lon, double h, long dat, double* field ) r = sqrt(r); - /* r is geocentric radial distance */ + /* r is geocentric radial distance */ c = cos(theta); s = sin(theta); - /*zero out arrays */ + /*zero out arrays */ for ( n = 0; n <= nmax; n++ ) { for ( m = 0; m <= n; m++ ) { P[n][m] = 0; @@ -271,7 +429,7 @@ double SGMagVar( double lat, double lon, double h, long dat, double* field ) field[4]=Y; field[5]=Z; /* output fields */ - /* find variation, leave in radians! */ + /* find variation, leave in radians! */ return atan2(Y, X); /* E is positive */ } - +#endif // TEST_NHV_HACKS