TTherminator/Therminator/Hypersurface.cxx

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