TTherminator/Therminator/Hypersurface.cxx

   1 #include <cstdlib>
   2 #include <cstring>
   3 #include "Hypersurface.h"
   4
   5 Hypersurface::Hypersurface(const char *dirname) {
   6 // Read FOHSI.txt file and init parameters
   7   char *dirfull;
   8   dirfull = (char *) malloc(sizeof(char)*strlen(dirname) + 5);
   9   sprintf(dirfull, "%s", dirname);
  10   if (strlen(dirfull) > 0)
  11     if (dirfull[strlen(dirfull)-1] != '/')
  12       strcat(dirfull, "/");
  13   FName = (char *) malloc(sizeof(char)* strlen(dirfull)+ 50);
  14   sprintf(FName,"%sFOHSI.txt", dirfull);
  15   HSFile = new ifstream;
  16   HSFile->open(FName);
  17   if(HSFile->is_open()) {
  18     HSFile->seekg (0, ios::beg);
  19 // phi
  20     HSFile->getline(buff,100);
  21     Np=atoi(buff);
  22     HSFile->getline(buff,100);
  23     ip=atof(buff);
  24     HSFile->getline(buff,100);
  25     fp=atof(buff);
  26 // zeta
  27     HSFile->getline(buff,100);
  28     Nz=atoi(buff);
  29     HSFile->getline(buff,100);
  30     iz=atof(buff);
  31     HSFile->getline(buff,100);
  32     fz=atof(buff);
  33 // freeze-out temperature
  34     HSFile->getline(buff,100);
  35     TFO=atof(buff);
  36 // hydro initial time
  37     HSFile->getline(buff,100);
  38     tau0=atof(buff);
  39
  40     HSFile->close();
  41     dp = (fp-ip)/(Np-1);
  42     dz = (fz-iz)/(Nz-1);
  43
  44 // create 2D array
  45     aArr   = new double* [Np];
  46     vArr   = new double* [Np];
  47     dArr   = new double* [Np];
  48     DpdArr = new double* [Np];
  49     DzdArr = new double* [Np];
  50     for(i=0;i<Np;i++) {
  51       aArr[i]   = new double [Nz];
  52       vArr[i]   = new double [Nz];
  53       dArr[i]   = new double [Nz];
  54       DpdArr[i] = new double [Nz];
  55       DzdArr[i] = new double [Nz];
  56     }
  57   } else {
  58     Np = 0; ip = 0.0; fp = 0.0; dp = 1.0;
  59     Nz = 0; iz = 0.0; fz = 0.0; dz = 1.0;
  60     TFO = 0.0;
  61     aArr   = NULL;
  62     vArr   = NULL;
  63     dArr   = NULL;
  64     DpdArr = NULL;
  65     DzdArr = NULL;
  66   }
  67   delete HSFile;
  68
  69   // Read FOHSa.txt file and fill array
  70   sprintf(FName,"%sFOHSa.txt", dirfull);
  71   //  FName = "FOHSa.txt";
  72   HSFile = new ifstream;
  73   HSFile->open(FName);
  74   if(HSFile->is_open()) {
  75     HSFile->seekg (0, ios::beg);
  76     for(i=0;i<Np;i++) {
  77       for(j=0;j<Nz;j++) {
  78         HSFile->getline(buff,100);
  79         aArr[i][j]=atof(buff);
  80       }
  81     }
  82     HSFile->close();
  83   }
  84   delete HSFile;
  85
  86   // Read FOHSv.txt file and fill array
  87   sprintf(FName,"%sFOHSv.txt", dirfull);
  88   //  FName = "FOHSv.txt";
  89   HSFile = new ifstream;
  90   HSFile->open(FName);
  91   if(HSFile->is_open()) {
  92     HSFile->seekg (0, ios::beg);
  93     for(i=0;i<Np;i++) {
  94       for(j=0;j<Nz;j++) {
  95         HSFile->getline(buff,100);
  96         vArr[i][j]=atof(buff);
  97       }
  98     }
  99     HSFile->close();
 100   }
 101   delete HSFile;
 102
 103   // Read FOHSd.txt file and fill array
 104   sprintf(FName,"%sFOHSd.txt", dirfull);
 105   //  FName = "FOHSd.txt";
 106   HSFile = new ifstream;
 107   HSFile->open(FName);
 108   if(HSFile->is_open()) {
 109     HSFile->seekg (0, ios::beg);
 110     for(i=0;i<Np;i++) {
 111       for(j=0;j<Nz;j++) {
 112         HSFile->getline(buff,100);
 113         dArr[i][j]=atof(buff);
 114       }
 115     }
 116     HSFile->close();
 117   }
 118   delete HSFile;
 119
 120   // Read FOHSDpd.txt file and fill array
 121   sprintf(FName,"%sFOHSDpd.txt", dirfull);
 122   //  FName = "FOHSDpd.txt";
 123   HSFile = new ifstream;
 124   HSFile->open(FName);
 125   if(HSFile->is_open()) {
 126     HSFile->seekg (0, ios::beg);
 127     for(i=0;i<Np;i++) {
 128       for(j=0;j<Nz;j++) {
 129         HSFile->getline(buff,100);
 130         DpdArr[i][j]=atof(buff);
 131       }
 132     }
 133     HSFile->close();
 134   }
 135   delete HSFile;
 136
 137   // Read FOHSDzd.txt file and fill array
 138   sprintf(FName,"%sFOHSDzd.txt", dirfull);
 139   //  FName = "FOHSDzd.txt";
 140   HSFile = new ifstream;
 141   HSFile->open(FName);
 142   if(HSFile->is_open()) {
 143     HSFile->seekg (0, ios::beg);
 144     for(i=0;i<Np;i++) {
 145       for(j=0;j<Nz;j++) {
 146         HSFile->getline(buff,100);
 147         DzdArr[i][j]=atof(buff);
 148       }
 149     }
 150     HSFile->close();
 151   }
 152   delete HSFile;
 153   free (FName);
 154   free (dirfull);
 155 }
 156
 157 Hypersurface::Hypersurface(void) {
 158   Hypersurface("./");
 159 }
 160
 161 Hypersurface::Hypersurface(const Hypersurface &aSurf)
 162 {
 163   i = aSurf.i;
 164   Hypersurface("./");
 165
 166 }
 167
 168 Hypersurface& Hypersurface::operator=(const Hypersurface &aSurf)
 169 {
 170   if (this != &aSurf) {
 171     i = aSurf.i;
 172     Hypersurface("./");
 173   }
 174
 175   return *this;
 176 }
 177
 178 Hypersurface::~Hypersurface(void) {
 179   for(i=0;i<Np;i++) {
 180     delete[] aArr[i];
 181     delete[] vArr[i];
 182     delete[] dArr[i];
 183     delete[] DpdArr[i];
 184     delete[] DzdArr[i];
 185   }
 186   delete[] aArr;
 187   delete[] vArr;
 188   delete[] dArr;
 189   delete[] DpdArr;
 190   delete[] DzdArr;
 191 }
 192
 193 double Hypersurface::fahs(double p, double z) {
 194   i=int(p/dp);
 195   j=int(z/dz);
 196   if(i>=Np-1) i=Np-2;
 197   if(j>=Nz-1) j=Nz-2;
 198   return (aArr[i][j]   * (i+1-p/dp) + aArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 199          (aArr[i][j+1] * (i+1-p/dp) + aArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 200 }
 201
 202 double Hypersurface::fvhs(double p, double z) {
 203   i=int(p/dp);
 204   j=int(z/dz);
 205   if(i>=Np-1) i=Np-2;
 206   if(j>=Nz-1) j=Nz-2;
 207   return (vArr[i][j]   * (i+1-p/dp) + vArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 208          (vArr[i][j+1] * (i+1-p/dp) + vArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 209 }
 210
 211 double Hypersurface::fdhs(double p, double z) {
 212   i=int(p/dp);
 213   j=int(z/dz);
 214   if(i>=Np-1) i=Np-2;
 215   if(j>=Nz-1) j=Nz-2;
 216   return (dArr[i][j]   * (i+1-p/dp) + dArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 217          (dArr[i][j+1] * (i+1-p/dp) + dArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 218 }
 219
 220 double Hypersurface::fDpdhs(double p, double z) {
 221   i=int(p/dp);
 222   j=int(z/dz);
 223   if(i>=Np-1) i=Np-2;
 224   if(j>=Nz-1) j=Nz-2;
 225   return (DpdArr[i][j]   * (i+1-p/dp) + DpdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 226          (DpdArr[i][j+1] * (i+1-p/dp) + DpdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 227 }
 228
 229 double Hypersurface::fDzdhs(double p, double z) {
 230   i=int(p/dp);
 231   j=int(z/dz);
 232   if(i>=Np-1) i=Np-2;
 233   if(j>=Nz-1) j=Nz-2;
 234   return (DzdArr[i][j]   * (i+1-p/dp) + DzdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 235          (DzdArr[i][j+1] * (i+1-p/dp) + DzdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 236 }