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::Hypersurface("./");
 159 }
 160
 161
 162 Hypersurface::~Hypersurface(void) {
 163   for(i=0;i<Np;i++) {
 164     delete[] aArr[i];
 165     delete[] vArr[i];
 166     delete[] dArr[i];
 167     delete[] DpdArr[i];
 168     delete[] DzdArr[i];
 169   }
 170   delete[] aArr;
 171   delete[] vArr;
 172   delete[] dArr;
 173   delete[] DpdArr;
 174   delete[] DzdArr;
 175 }
 176
 177 double Hypersurface::fahs(double p, double z) {
 178   i=int(p/dp);
 179   j=int(z/dz);
 180   if(i>=Np-1) i=Np-2;
 181   if(j>=Nz-1) j=Nz-2;
 182   return (aArr[i][j]   * (i+1-p/dp) + aArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 183          (aArr[i][j+1] * (i+1-p/dp) + aArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 184 }
 185
 186 double Hypersurface::fvhs(double p, double z) {
 187   i=int(p/dp);
 188   j=int(z/dz);
 189   if(i>=Np-1) i=Np-2;
 190   if(j>=Nz-1) j=Nz-2;
 191   return (vArr[i][j]   * (i+1-p/dp) + vArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 192          (vArr[i][j+1] * (i+1-p/dp) + vArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 193 }
 194
 195 double Hypersurface::fdhs(double p, double z) {
 196   i=int(p/dp);
 197   j=int(z/dz);
 198   if(i>=Np-1) i=Np-2;
 199   if(j>=Nz-1) j=Nz-2;
 200   return (dArr[i][j]   * (i+1-p/dp) + dArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 201          (dArr[i][j+1] * (i+1-p/dp) + dArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 202 }
 203
 204 double Hypersurface::fDpdhs(double p, double z) {
 205   i=int(p/dp);
 206   j=int(z/dz);
 207   if(i>=Np-1) i=Np-2;
 208   if(j>=Nz-1) j=Nz-2;
 209   return (DpdArr[i][j]   * (i+1-p/dp) + DpdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 210          (DpdArr[i][j+1] * (i+1-p/dp) + DpdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 211 }
 212
 213 double Hypersurface::fDzdhs(double p, double z) {
 214   i=int(p/dp);
 215   j=int(z/dz);
 216   if(i>=Np-1) i=Np-2;
 217   if(j>=Nz-1) j=Nz-2;
 218   return (DzdArr[i][j]   * (i+1-p/dp) + DzdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
 219          (DzdArr[i][j+1] * (i+1-p/dp) + DzdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
 220 }