vhlle/src/s95p.cpp

0001 #include <cmath>
0002 #include <iostream>
0003 #include <fstream>
0004 #include <cstdlib>
0005
0006 extern int glauberVariable;
0007
0008 namespace s95p {
0009 using namespace std;
0010
0011 int NX, NY;  // N*N table
0012 double** e, xmin, xmax, ymin, ymax;
0013
0014 double s95p_s(double e) {
0015   double val;
0016   if (e < 0.1270769021427449)
0017     val = 12.2304 * pow(e, 1.16849);
0018   else if (e < 0.4467079524674040)
0019     val = 11.9279 * pow(e, 1.15635);
0020   else if (e < 1.9402832534193788)
0021     val = 0.0580578 + 11.833 * pow(e, 1.16187);
0022   else if (e < 3.7292474570977285)
0023     val =
0024         18.202 * e - 62.021814 - 4.85479 * exp(-2.72407e-11 * pow(e, 4.54886)) +
0025         65.1272 * pow(e, -0.128012) * exp(-0.00369624 * pow(e, 1.18735)) -
0026         4.75253 * pow(e, -1.18423);
0027   else
0028     val = 18.202 * e - 63.0218 - 4.85479 * exp(-2.72407e-11 * pow(e, 4.54886)) +
0029           65.1272 * pow(e, -0.128012) * exp(-0.00369624 * pow(e, 1.18735));
0030   return pow(val, 3. / 4.);
0031 }
0032
0033 double s95p_e(double s) {
0034   double e = pow(s, 4. / 3.) / 18.2, e0;
0035   do {
0036     e0 = e;
0037     e = e0 - (s95p_s(e0) - s) / (6.608681233 * pow(e0, -0.25));
0038   } while (fabs(e - e0) > 1e-10);
0039   return e;
0040 }
0041
0042 void loadSongIC(char* filename, double factor) {
0043   ifstream fin(filename);
0044   if (!fin) {
0045     cout << "cannot open " << filename << endl;
0046     exit(1);
0047   }
0048   fin >> NX >> NY;
0049   double* x = new double[NX];
0050   double* y = new double[NY];
0051   e = new double* [NX];
0052   for (int ix = 0; ix < NX; ix++) e[ix] = new double[NY];
0053
0054   for (int ix = 0; ix < NX; ix++)
0055     for (int iy = 0; iy < NY; iy++) {
0056       fin >> x[ix] >> y[iy] >> e[ix][iy];
0057       if (glauberVariable == 1)
0058         e[ix][iy] = s95p_e(factor * e[ix][iy]);
0059       else
0060         e[ix][iy] = factor * e[ix][iy];
0061     }
0062   xmin = x[0];
0063   xmax = x[NX - 1];
0064   ymin = y[0];
0065   ymax = y[NY - 1];
0066   delete[] x;
0067   delete[] y;
0068 }
0069
0070 double getSongEps(double x, double y) {
0071   // linear interpolation in 2D is used
0072   const double dx = (xmax - xmin) / (NX - 1);
0073   const double dy = (ymax - ymin) / (NY - 1);
0074   const int ix = (int)((x - xmin) / dx);
0075   const int iy = (int)((y - ymin) / dy);
0076   if (ix < 0 || ix > NX - 2 || iy < 0 || iy > NY - 2) return 0.0;
0077   const double sx = x - xmin - ix * dx;
0078   const double sy = y - ymin - iy * dy;
0079   double wx[2] = {(1. - sx / dx), sx / dx};
0080   double wy[2] = {(1. - sy / dy), sy / dy};
0081   double result = 0;
0082   for (int jx = 0; jx < 2; jx++)
0083     for (int jy = 0; jy < 2; jy++) {
0084       result += wx[jx] * wy[jy] * e[ix + jx][iy + jy];
0085     }
0086   return result;
0087 }
0088
0089 }  // end s95p