update Voigt2d ABC regions for comparison

define.h

@@ -36,7 +36,7 @@
 #define def_T0 296.0 		//Reference Temperature in K
 #define def_PObar 0.986923	//Referecne Pressure 1 bar in atm for ExoMol
 #define def_POatm 1.0		//Referecne Pressure 1 atm in atm for Hitran
-#define def_kB 1.3806489e-16 	//Boltzmann constant in erg/K
+#define def_kB 1.3806489e-16 	//Boltzmann constant in erg/K = cm^2 g / (s^2 K)
 #define def_h 6.62606957e-27	//Planck costant in erg s
 #define def_c 2.99792458e10 	//Speed of light cm/s
 #define def_NA 6.0221412927e23	//Avogadro Constant  1/mol

heliosk.cu

@@ -179,8 +179,8 @@ int main(int argc, char*argv[]){
 double xMax = 10.0;
 double yMax = 10.0;
 
-int Nx = 10000;
-int Ny = 10000;
+int Nx = 1000;
+int Ny = 1000;
 
 int Nxtex = Nx + 1;
 int Nytex = Ny + 1;
@@ -200,20 +200,25 @@ cudaMallocPitch((void **) &K2d_d, &pitch, Nxtex * sizeof(double), Nytex);
 	double a = (double)(M_PI * sqrt(-1.0 / log(def_TOLF * 0.5)));
 	double b = (double)(1.0 / sqrt(M_PI));
 	double c = (double)(2.0 * a / M_PI);
+
 	Voigt_2d_kernel <<< dim3((Nxtex + 31) / 32, (Nytex + 31) / 32), dim3(32, 32, 1) >>> (a, b, c, K2d_d, Nxtex, Nytex, pitch, xMax, xMax);
 	cudaMemcpy2D(K2d_h, Nxtex * sizeof(double), K2d_d, pitch, Nxtex * sizeof(double), Nytex, cudaMemcpyDeviceToHost);
 }
-/ *
+// / *
 for(int i = 0; i < Nxtex - 1; ++i){
 	for(int j = 0; j < Nytex - 1; ++j){
-		double x = i * xMax / double(Nxtex);
-		double y = j * yMax / double(Nytex);
-		if( x < xMax && y < yMax){
+		//x and y arrays from 0.1 to 2000
+		double x = exp(-2.3 + i * xMax / double(Nxtex - 1));
+		double y = exp(-2.3 + j * yMax / double(Nytex - 1));
+		//if( x < xMax && y < yMax){
 			printf("%g %g %.15g\n", x, y, K2d_h[j * Nxtex + i]);
-		}
+		//}
 	}
 }
-* /
+// * /
+return 0;
+}
+/ *
 
 float *K2df_h, *K2df_d;
 size_t pitchf;

prepareExomol.cpp

@@ -210,8 +210,8 @@ int main(int argc, char*argv[]){
 	char qFilename[15][160];
 	Init(m, param, qFilename);
 
-	double mass = m.ISO[0].m;	//Molar Mass (g)
-	mass /= def_NA;
+	double mass = m.ISO[0].m;	//Molar Mass (g/mol)
+	mass /= def_NA;			//mass in g
 	int *id, *g;
 	double *E;
 

voigt.h

@@ -568,8 +568,9 @@ __global__ void Voigt_2d_kernel(const double a, const double b, const double c,
 	int idy = blockIdx.y * blockDim.y + threadIdx.y;
 
 	if(idx < Nx && idy < Ny){
-		double x = idx * xMax / double(Nx - 1);
-		double y = idy * yMax / double(Ny - 1);
+		//x and y arrays from 0.1 to 2000
+		double x = exp(-2.3 + idx * xMax / double(Nx - 1));
+		double y = exp(-2.3 + idy * yMax / double(Ny - 1));
 		double s1, s2, s3;
 		double ex2 = expf(-x * x);
 
@@ -589,10 +590,25 @@ __global__ void Voigt_2d_kernel(const double a, const double b, const double c,
 		if(y == 0.0) t1 = ex2;
 
 		//K_d[idy * Nx + idx] = t1 * b;
+
+		double xxyy = x * x + y * y;
 		double *row = (double *)(((char *)K_d)+(idy*pitch));
-		row[idx] = t1 * b;
-//if(idy == 0) printf("a %d %d %g %g %g\n", idx, idy, x, y, float(idy * Nx + idx));
-//printf("%g %g %g %g %g %g\n", x, y, s1, s2, s3, K_d[idy * Nx + idx]);
+		if(xxyy < 100.0){
+			row[idx] = t1 * b;
+		}
+		else if(xxyy < 1.0e6){
+			//Region B
+			double t1 = x * x * 6.0;
+			double t2 = xxyy + 1.5;
+			double t3 = M_PI * t2; 
+			double t4 = (t3 * (2.0 * t2 + xxyy) - 2.0 * t1) / (3.0 * xxyy * (t3 * t2 - t1));
+			row[idx] = y * t4 / M_PI;
+
+		}
+		else{
+			//Region A
+			row[idx] = y / (M_PI * xxyy);
+		}
 	}
 }
 __global__ void Voigt_2df_kernel(const float a, const float b, const float c, float *K_d, int Nx, int Ny, size_t pitch, float xMax, float yMax){