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