Albany: a Trilinos-based PDE code: AAdapt_AnalyticFunction.cpp Source File

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <cmath>
00008 #include <ctime>
00009 #include <cstdlib>
00010 
00011 #include "AAdapt_AnalyticFunction.hpp"
00012 #include "Teuchos_TestForException.hpp"
00013 #include "Teuchos_Exceptions.hpp"
00014 
00015 const double pi = 3.141592653589793;
00016 
00017 
00018 // Factory method to build functions based on a string
00019 Teuchos::RCP<AAdapt::AnalyticFunction> AAdapt::createAnalyticFunction(
00020   std::string name, int neq, int numDim,
00021   Teuchos::Array<double> data) {
00022   Teuchos::RCP<AAdapt::AnalyticFunction> F;
00023 
00024   if(name == "Constant")
00025     F = Teuchos::rcp(new AAdapt::ConstantFunction(neq, numDim, data));
00026 
00027   else if(name == "1D Gauss-Sin")
00028     F = Teuchos::rcp(new AAdapt::GaussSin(neq, numDim, data));
00029 
00030   else if(name == "1D Gauss-Cos")
00031     F = Teuchos::rcp(new AAdapt::GaussCos(neq, numDim, data));
00032 
00033   else if(name == "Linear Y")
00034     F = Teuchos::rcp(new AAdapt::LinearY(neq, numDim, data));
00035 
00036   else if(name == "Gaussian Pressure")
00037     F = Teuchos::rcp(new AAdapt::GaussianPress(neq, numDim, data));
00038 
00039   else if(name == "Sin-Cos")
00040     F = Teuchos::rcp(new AAdapt::SinCos(neq, numDim, data));
00041 
00042   else if(name == "Taylor-Green Vortex")
00043     F = Teuchos::rcp(new AAdapt::TaylorGreenVortex(neq, numDim, data));
00044 
00045   else if(name == "1D Acoustic Wave")
00046     F = Teuchos::rcp(new AAdapt::AcousticWave(neq, numDim, data));
00047 
00048   else if(name == "Aeras Schar Density")
00049     F = Teuchos::rcp(new AAdapt::AerasScharDensity(neq, numDim, data));
00050 
00051   else if(name == "Aeras Heaviside")
00052     F = Teuchos::rcp(new AAdapt::AerasHeaviside(neq, numDim, data));
00053 
00054   else if(name == "Aeras CosineBell")
00055       F = Teuchos::rcp(new AAdapt::AerasCosineBell(neq, numDim, data));
00056 
00057   else if(name == "Aeras ZonalFlow") //this used to be called TestCase2.  Irina has renamed it so it can be used for test case 5 too.
00058       F = Teuchos::rcp(new AAdapt::AerasZonalFlow(neq, numDim, data));
00059 
00060   else if(name == "Aeras PlanarCosineBell")
00061         F = Teuchos::rcp(new AAdapt::AerasPlanarCosineBell(neq, numDim, data));
00062 
00063   else if(name == "Aeras RossbyHaurwitzWave")
00064       F = Teuchos::rcp(new AAdapt::AerasRossbyHaurwitzWave(neq, numDim, data));
00065 
00066    else
00067     TEUCHOS_TEST_FOR_EXCEPTION(name != "Valid Initial Condition Function",
00068                                std::logic_error,
00069                                "Unrecognized initial condition function name: " << name);
00070 
00071   return F;
00072 }
00073 
00074 
00075 //*****************************************************************************
00076 AAdapt::ConstantFunction::ConstantFunction(int neq_, int numDim_,
00077     Teuchos::Array<double> data_) : numDim(numDim_), neq(neq_), data(data_) {
00078   TEUCHOS_TEST_FOR_EXCEPTION((data.size() != neq),
00079                              std::logic_error,
00080                              "Error! Invalid specification of initial condition: incorrect length of Function Data for Constant Function; neq = " << neq << ", data.size() = " << data.size() <<  std::endl) ;
00081 }
00082 void AAdapt::ConstantFunction::compute(double* x, const double* X) {
00083   if(data.size() > 0)
00084     for(int i = 0; i < neq; i++)
00085       x[i] = data[i];
00086 }
00087 
00088 //*****************************************************************************
00089 // Private convenience function
00090 long AAdapt::seedgen(int worksetID) {
00091   long seconds, s, seed, pid;
00092 
00093   pid = getpid();
00094   s = time(&seconds);    /* get CPU seconds since 01/01/1970 */
00095 
00096   // Use worksetID to give more randomness between calls
00097 
00098   seed = abs(((s * 181) * ((pid - 83) * 359) * worksetID) % 104729);
00099   return seed;
00100 }
00101 
00102 AAdapt::ConstantFunctionPerturbed::ConstantFunctionPerturbed(int neq_, int numDim_,
00103     int worksetID,
00104     Teuchos::Array<double> data_,  Teuchos::Array<double> pert_mag_)
00105   : numDim(numDim_), neq(neq_), data(data_), pert_mag(pert_mag_) {
00106 
00107   TEUCHOS_TEST_FOR_EXCEPTION((data.size() != neq || pert_mag.size() != neq),
00108                              std::logic_error,
00109                              "Error! Invalid specification of initial condition: incorrect length of " <<
00110                              "Function Data for Constant Function Perturbed; neq = " << neq <<
00111                              ", data.size() = " << data.size()
00112                              << ", pert_mag.size() = " << pert_mag.size()
00113                              <<  std::endl) ;
00114 
00115   //  srand( time(NULL) ); // seed the random number gen
00116   srand(seedgen(worksetID)); // seed the random number gen
00117 
00118 }
00119 
00120 void AAdapt::ConstantFunctionPerturbed::compute(double* x, const double* X) {
00121   for(int i = 0; i < neq; i++)
00122     x[i] = data[i] + udrand(-pert_mag[i], pert_mag[i]);
00123 }
00124 
00125 // Private convenience function
00126 double AAdapt::ConstantFunctionPerturbed::udrand(double lo, double hi) {
00127   static const double base = 1.0 / (RAND_MAX + 1.0);
00128   double deviate = std::rand() * base;
00129   return lo + deviate * (hi - lo);
00130 }
00131 //*****************************************************************************
00132 AAdapt::ConstantFunctionGaussianPerturbed::ConstantFunctionGaussianPerturbed(int neq_, int numDim_,
00133     int worksetID,
00134     Teuchos::Array<double> data_,  Teuchos::Array<double> pert_mag_)
00135   : numDim(numDim_),
00136     neq(neq_),
00137     data(data_),
00138     pert_mag(pert_mag_),
00139     //      rng(boost::random::random_device()()), // seed the rng
00140     rng(seedgen(worksetID)), // seed the rng
00141     nd(neq_),
00142     var_nor(neq_) {
00143 
00144 
00145   TEUCHOS_TEST_FOR_EXCEPTION((data.size() != neq || pert_mag.size() != neq),
00146                              std::logic_error,
00147                              "Error! Invalid specification of initial condition: incorrect length of " <<
00148                              "Function Data for Constant Function Gaussian Perturbed; neq = " << neq <<
00149                              ", data.size() = " << data.size()
00150                              << ", pert_mag.size() = " << pert_mag.size()
00151                              <<  std::endl) ;
00152 
00153   if(data.size() > 0 && pert_mag.size() > 0)
00154     for(int i = 0; i < neq; i++)
00155       if(pert_mag[i] > std::numeric_limits<double>::epsilon()) {
00156 
00157         nd[i] = Teuchos::rcp(new boost::normal_distribution<double>(data[i], pert_mag[i]));
00158         var_nor[i] = Teuchos::rcp(new
00159                                   boost::variate_generator<boost::mt19937&, boost::normal_distribution<double> >(rng, *nd[i]));
00160 
00161       }
00162 
00163 }
00164 
00165 void AAdapt::ConstantFunctionGaussianPerturbed::compute(double* x, const double* X) {
00166 
00167   for(int i = 0; i < neq; i++)
00168     if(var_nor[i] != Teuchos::null)
00169       x[i] = (*var_nor[i])();
00170 
00171     else
00172       x[i] = data[i];
00173 
00174 }
00175 
00176 
00177 //*****************************************************************************
00178 AAdapt::GaussSin::GaussSin(int neq_, int numDim_, Teuchos::Array<double> data_)
00179   : numDim(numDim_), neq(neq_), data(data_) {
00180   TEUCHOS_TEST_FOR_EXCEPTION((neq != 1) || (numDim != 1) || (data.size() != 1),
00181                              std::logic_error,
00182                              "Error! Invalid call of GaussSin with " << neq
00183                              << " " << numDim << "  " << data.size() << std::endl);
00184 }
00185 void AAdapt::GaussSin::compute(double* x, const double* X) {
00186   x[0] =     sin(pi * X[0]) + 0.5 * data[0] * X[0] * (1.0 - X[0]);
00187 }
00188 
00189 //*****************************************************************************
00190 AAdapt::GaussCos::GaussCos(int neq_, int numDim_, Teuchos::Array<double> data_)
00191   : numDim(numDim_), neq(neq_), data(data_) {
00192   TEUCHOS_TEST_FOR_EXCEPTION((neq != 1) || (numDim != 1) || (data.size() != 1),
00193                              std::logic_error,
00194                              "Error! Invalid call of GaussCos with " << neq
00195                              << " " << numDim << "  " << data.size() << std::endl);
00196 }
00197 void AAdapt::GaussCos::compute(double* x, const double* X) {
00198   x[0] = 1 + cos(2 * pi * X[0]) + 0.5 * data[0] * X[0] * (1.0 - X[0]);
00199 }
00200 //*****************************************************************************
00201 AAdapt::LinearY::LinearY(int neq_, int numDim_, Teuchos::Array<double> data_)
00202   : numDim(numDim_), neq(neq_), data(data_) {
00203   TEUCHOS_TEST_FOR_EXCEPTION((neq < 2) || (numDim < 2) || (data.size() != 1),
00204                              std::logic_error,
00205                              "Error! Invalid call of LinearY with " << neq
00206                              << " " << numDim << "  " << data.size() << std::endl);
00207 }
00208 void AAdapt::LinearY::compute(double* x, const double* X) {
00209   x[0] = 0.0;
00210   x[1] = data[0] * X[0];
00211 
00212   if(numDim > 2) x[2] = 0.0;
00213 }
00214 //*****************************************************************************
00215 AAdapt::GaussianPress::GaussianPress(int neq_, int numDim_, Teuchos::Array<double> data_)
00216   : numDim(numDim_), neq(neq_), data(data_) {
00217   TEUCHOS_TEST_FOR_EXCEPTION((neq < 3) || (numDim < 2) || (data.size() != 4),
00218                              std::logic_error,
00219                              "Error! Invalid call of GaussianPress with " << neq
00220                              << " " << numDim << "  " << data.size() << std::endl);
00221 }
00222 void AAdapt::GaussianPress::compute(double* x, const double* X) {
00223   for(int i = 0; i < neq - 1; i++) {
00224     x[i] = 0.0;
00225   }
00226 
00227   x[neq - 1] = data[0] * exp(-data[1] * ((X[0] - data[2]) * (X[0] - data[2]) + (X[1] - data[3]) * (X[1] - data[3])));
00228 }
00229 //*****************************************************************************
00230 AAdapt::SinCos::SinCos(int neq_, int numDim_, Teuchos::Array<double> data_)
00231   : numDim(numDim_), neq(neq_), data(data_) {
00232   TEUCHOS_TEST_FOR_EXCEPTION((neq < 3) || (numDim < 2),
00233                              std::logic_error,
00234                              "Error! Invalid call of SinCos with " << neq
00235                              << " " << numDim << "  " << data.size() << std::endl);
00236 }
00237 void AAdapt::SinCos::compute(double* x, const double* X) {
00238   x[0] = sin(2.0 * pi * X[0]) * cos(2.0 * pi * X[1]);
00239   x[1] = cos(2.0 * pi * X[0]) * sin(2.0 * pi * X[1]);
00240   x[2] = sin(2.0 * pi * X[0]) * sin(2.0 * pi * X[1]);
00241 }
00242 //*****************************************************************************
00243 AAdapt::TaylorGreenVortex::TaylorGreenVortex(int neq_, int numDim_, Teuchos::Array<double> data_)
00244   : numDim(numDim_), neq(neq_), data(data_) {
00245   TEUCHOS_TEST_FOR_EXCEPTION((neq < 3) || (numDim != 2),
00246                              std::logic_error,
00247                              "Error! Invalid call of TaylorGreenVortex with " << neq
00248                              << " " << numDim << "  " << data.size() << std::endl);
00249 }
00250 void AAdapt::TaylorGreenVortex::compute(double* x, const double* X) {
00251   x[0] = 1.0; //initial density
00252   x[1] = -cos(2.0 * pi * X[0]) * sin(2.0 * pi * X[1]); //initial u-velocity
00253   x[2] = sin(2.0 * pi * X[0]) * cos(2.0 * pi * X[1]); //initial v-velocity
00254   x[3] = cos(2.0 * pi * X[0]) + cos(2.0 * pi * X[1]); //initial temperature
00255 }
00256 //*****************************************************************************
00257 AAdapt::AcousticWave::AcousticWave(int neq_, int numDim_, Teuchos::Array<double> data_)
00258   : numDim(numDim_), neq(neq_), data(data_) {
00259   TEUCHOS_TEST_FOR_EXCEPTION((neq > 3) || (numDim > 2) || (data.size() != 3),
00260                              std::logic_error,
00261                              "Error! Invalid call of AcousticWave with " << neq
00262                              << " " << numDim << "  " << data.size() << std::endl);
00263 }
00264 void AAdapt::AcousticWave::compute(double* x, const double* X) {
00265   const double U0 = data[0];
00266   const double n = data[1];
00267   const double L = data[2];
00268   x[0] = U0 * cos(n * X[0] / L);
00269 
00270   for(int i = 1; i < numDim; i++)
00271     x[i] = 0.0;
00272 }
00273 
00274 //*****************************************************************************
00275 AAdapt::AerasScharDensity::AerasScharDensity(int neq_, int numDim_, Teuchos::Array<double> data_)
00276   : numDim(numDim_), neq(neq_), data(data_) {
00277   TEUCHOS_TEST_FOR_EXCEPTION((neq > 1) || (numDim > 2),
00278                              std::logic_error,
00279                              "Error! Invalid call of Aeras Schar Density with " << neq
00280                              << " " << numDim << std::endl);
00281 }
00282 void AAdapt::AerasScharDensity::compute(double* x, const double* X) {
00283   //const double U0 = data[0];
00284   double r = sqrt ( std::pow((X[0] - 100.0)/25.0 ,2) +  std::pow((X[1] - 9.0)/3.0,2));
00285   if (r <= 1.0) x[0] = std::pow(cos(pi*r / 2.0),2);
00286   else          x[0] = 0.0;
00287 }
00288 //*****************************************************************************
00289 AAdapt::AerasHeaviside::AerasHeaviside(int neq_, int numDim_, Teuchos::Array<double> data_)
00290   : numDim(numDim_), neq(neq_), data(data_) {
00291   TEUCHOS_TEST_FOR_EXCEPTION((neq != 3) || (numDim != 2),
00292                              std::logic_error,
00293                              "Error! Invalid call of Aeras Heaviside with " << neq
00294                              << " " << numDim <<  std::endl);
00295 }
00296 void AAdapt::AerasHeaviside::compute(double* x, const double* X) {
00297   //const double U0 = data[0];
00298   if (X[0] <= 0.5) x[0] = 1.1;
00299   else             x[0] = 1.0;
00300   x[1]=0.0;
00301   x[2]=0.0;
00302 }
00303 //*****************************************************************************
00304 AAdapt::AerasCosineBell::AerasCosineBell(int neq_, int spatialDim_, Teuchos::Array<double> data_)
00305   : spatialDim(spatialDim_), neq(neq_), data(data_) {
00306   TEUCHOS_TEST_FOR_EXCEPTION( (neq!=3 || spatialDim!=3 || data.size()!=2) ,
00307                              std::logic_error,
00308                              "Error! Invalid call of Aeras CosineBell with " << neq
00309                              << " " << spatialDim <<  " "<< data.size()<< std::endl);
00310 }
00311 void AAdapt::AerasCosineBell::compute(double* solution, const double* X) {
00312   const double u0 = data[0];  // magnitude of wind
00313   const double cosAlpha = std::cos(data[1]);
00314   const double sinAlpha = std::sin(data[1]);
00315 
00316   const double lambda_c = 1.5*pi;
00317   const double theta_c = 0;
00318   const double sinTheta_c = std::sin(theta_c);
00319   const double cosTheta_c = std::cos(theta_c);
00320 
00321   const double x = X[0];
00322   const double y = X[1];
00323   const double z = X[2];
00324 
00325   const double sinTheta = z;
00326 //  const double cosTheta = std::sqrt(x*x + y*y);
00327 
00328   const double cosTheta = std::sqrt(1-sinTheta*sinTheta);
00329   const double sinLambda = std::abs(cosTheta)>.0001 ? y/cosTheta : 0;
00330   const double cosLambda = std::abs(cosTheta)>.0001 ? x/cosTheta : 1;
00331 
00332   const double u = u0*(cosTheta*cosAlpha + sinTheta*cosLambda*sinAlpha);
00333   const double v = -u0*(sinLambda*sinAlpha);
00334 
00335   const double a = 1; //radius of earth;
00336   const double R = a/3.;
00337   const double h0 = 1000./6378100.0;   // 1000/radius o earth in meters
00338 
00339   //const double lambda =  std::atan2(y,x);
00340   const double theta  = std::asin(sinTheta);
00341   static const double DIST_THRESHOLD = 1.0e-9;
00342   double lambda = std::atan2(y, x);
00343   if (std::abs(std::abs(theta)-pi/2) < DIST_THRESHOLD) lambda = 0;
00344   else if (lambda < 0) lambda += 2*pi;
00345 
00346   const double r = a*std::acos(sinTheta_c*sinTheta + cosTheta_c*cosTheta*std::cos(lambda - lambda_c));
00347 
00348   const double h = r < R ? 0.5*h0*(1 + std::cos(pi*r/R)) : 0;
00349 
00350   solution[0] = h;
00351   solution[1] = u;
00352   solution[2] = v;
00353 }
00354 //*****************************************************************************
00355 //IK, 2/5/14: added to data array h0*g, which corresponds to data[2] 
00356 AAdapt::AerasZonalFlow::AerasZonalFlow(int neq_, int spatialDim_, Teuchos::Array<double> data_)
00357   : spatialDim(spatialDim_), neq(neq_), data(data_) {
00358   TEUCHOS_TEST_FOR_EXCEPTION( (neq!=3 || spatialDim!=3 || data.size()!=3) ,
00359                              std::logic_error,
00360                              "Error! Invalid call of Aeras ZonalFlow with " << neq
00361                              << " " << spatialDim <<  " "<< data.size()<< std::endl);
00362 }
00363 void AAdapt::AerasZonalFlow::compute(double* solution, const double* X) {
00364   const double u0 = data[0];  // magnitude of wind
00365   const double cosAlpha = std::cos(data[1]);  //alpha
00366   const double sinAlpha = std::sin(data[1]);
00367   const double h0g = data[3]; //h0*g 
00368 
00369   const double x = X[0];
00370   const double y = X[1];
00371   const double z = X[2];
00372 
00373   // ==========================================================
00374   // enforce three facts:
00375   //
00376   // 1) lon at poles is defined to be zero
00377   //
00378   // 2) Grid points must be separated by about .01 Meter (on earth)
00379   //   from pole to be considered "not the pole".
00380   //
00381   // 3) range of lon is { 0<= lon < 2*PI }
00382   //
00383   // ==========================================================
00384   static const double pi = 3.1415926535897932385;
00385   static const double DIST_THRESHOLD = 1.0e-9;
00386   const double Theta  = std::asin(z);
00387 
00388   double Lambda = std::atan2(y,x);
00389   if (std::abs(std::abs(Theta)-pi/2) < DIST_THRESHOLD) Lambda = 0;
00390   else if (Lambda < 0) Lambda += 2*pi;
00391 
00392   const double sinTheta  = std::sin(Theta);
00393   const double cosTheta  = std::cos(Theta);
00394   const double sinLambda = std::sin(Lambda);
00395   const double cosLambda = std::cos(Lambda);
00396 
00397   const double u = u0*(cosTheta*cosAlpha + sinTheta*cosLambda*sinAlpha);
00398   const double v = -u0*(sinLambda*sinAlpha);
00399 
00400   const double a      = 6.37122e06; //radius of earth (m);
00401   const double g      = 9.80616;    // gravity m/s/s
00402   const double h0     = h0g/g;  // 1000/radius o earth in (m)
00403   const double Omega  = 7.292e-05; // 1/sec.
00404 
00405   const double lambda =  std::atan2(x,y); //longitude 
00406 
00407   //IK, 2/5/14: it is best not to use pow, I think. 
00408   const double h = h0 - 1.0/g * (a*Omega*u0 + u0*u0/2.0)*(-1.0*cosLambda*cosTheta*sinAlpha + sinTheta*cosAlpha)*(-1.0*cosLambda*cosTheta*sinAlpha + sinTheta*cosAlpha);
00409 
00410   solution[0] = h;
00411   solution[1] = u;
00412   solution[2] = v;
00413 }
00414 //*****************************************************************************
00415 AAdapt::AerasPlanarCosineBell::AerasPlanarCosineBell(int neq_, int numDim_, Teuchos::Array<double> data_)
00416   : numDim(numDim_), neq(neq_), data(data_) {
00417   TEUCHOS_TEST_FOR_EXCEPTION( (neq!=3 || numDim!=2 || data.size()!=3)  ,
00418                              std::logic_error,
00419                              "Error! Invalid call of Aeras PlanarCosineBell with " << neq
00420                              << " " << numDim <<  " "<< data.size()<< std::endl);
00421 }
00422 void AAdapt::AerasPlanarCosineBell::compute(double* solution, const double* X) {
00423   const double u0 = data[0];  // magnitude of wind
00424   const double h0 = data[1];
00425   const double R = data[2];
00426 
00427 
00428   const double x = X[0];
00429   const double y = X[1];
00430   const double z = X[2];
00431 
00432   const double u = u0;
00433   const double v = 0;
00434 
00435 
00436   const double r = std::sqrt( (x-0.5)*(x-0.5) + (y-0.5)*(y-0.5));
00437 
00438   const double h = r < R ? 1 + 0.5*h0*(1 + std::cos(pi*r/R)) : 1;
00439 
00440   solution[0] = h;
00441   solution[1] = u;
00442   solution[2] = v;
00443 
00444 }
00445 //*****************************************************************************
00446 AAdapt::AerasRossbyHaurwitzWave::AerasRossbyHaurwitzWave(int neq_, int spatialDim_, Teuchos::Array<double> data_)
00447   : spatialDim(spatialDim_), neq(neq_), data(data_) {
00448   TEUCHOS_TEST_FOR_EXCEPTION( (neq!=3 || spatialDim!=3 || data.size()!=4) ,
00449                              std::logic_error,
00450                              "Error! Invalid call of Aeras RossbyHaurwitzWave with " << neq
00451                              << " " << spatialDim <<  " "<< data.size()<< std::endl);
00452 }
00453 void AAdapt::AerasRossbyHaurwitzWave::compute(double* solution, const double* X) {
00454 
00455   // Problem constants
00456   const double a     = 1;         // Radius of earth
00457   const double aSq   = a * a;     // Radius of earth squared
00458   const double g     = 9.80616;   // Gravitational acceleration
00459   const double Omega = 7.292e-5;  // Angular rotation of earth
00460 
00461   // User-supplied data
00462   const double omega = data[0];
00463   const double K     = data[1];
00464   const int    R     = data[2];
00465   const double h0    = data[3];
00466 
00467   // Computed constants
00468   const double KSq = K * K;
00469   const int    RSq = R * R;
00470 
00471   // Coordinates
00472   const double x = X[0];
00473   const double y = X[1];
00474   const double z = X[2];
00475 
00476   // Trigonometric forms of latitude and longitude
00477   const double lambda     = std::atan2(x,y);
00478   const double sinRLambda = std::sin(R*lambda);
00479   const double cosRLambda = std::cos(R*lambda);
00480   const double sinTheta   = z;
00481   const double cosTheta   = std::sqrt(x*x + y*y);
00482   const double cosSqTheta = cosTheta * cosTheta;
00483 
00484   // Initial velocities
00485   const double u = a * omega * cosTheta + a * K * std::pow(cosTheta,R-1) *
00486     (R * sinTheta*sinTheta - cosSqTheta) * cosRLambda;
00487   const double v = -a*K*R * std::pow(cosTheta,R-1) * sinTheta * sinRLambda;
00488 
00489   // Latitudinal constants for computing h
00490   const double A = 0.5 * omega * (2*Omega + omega) * cosSqTheta + 0.25 * KSq *
00491                    std::pow(cosSqTheta,R) + ((R+1) * cosSqTheta +
00492                    (2*RSq - R - 2) - 2 * RSq * std::pow(cosTheta,-2));
00493   const double B = (2 * (Omega + omega) * K * std::pow(cosTheta,R) *
00494                     ((RSq + 2*R + 2) - (R+1) * (R+1) * cosSqTheta)) /
00495                    ((R+1)*(R+2));
00496   const double C = 0.25 * KSq * std::pow(cosSqTheta,R) *
00497                    ((R+1) * cosSqTheta - (R+2));
00498 
00499   // Height field
00500   const double h = h0 +
00501                    (aSq*A + aSq*B*cosRLambda + aSq*C*std::cos(2*R*lambda)) / g;
00502 
00503   // Assign the solution
00504   solution[0] = h;
00505   solution[1] = u;
00506   solution[2] = v;
00507 }