Add option for external background point source

.gitignore

@@ -330,4 +330,7 @@ Codes/*.csv
 Codes/*.dat
 Codes/*.dem
 Codes/*.txt
-Codes/*.xlsx
+Codes/*.xlsx
+Codes/*.zip
+Codes/*.7z
+Codes/*.rar

Codes/coefficients_calculator.c

@@ -13,8 +13,10 @@
 // E[TRANS_N]           : Boltzmann factor, exp(-dEJJ'/kT), for energy difference between J and J', dEJJ'. EJJ' = E[(J-1)J/2+J']
 // F[LEVEL_N-1]         : Flux normalization factor, 2h(vJJ')^3/c^2, FJJ' = F[J']
 // Br_n[LEVEL_N-1]      : Normalized cosmic blackbody radiation intensity for J -> J'=J-1, Br_JJ' = Br[J']
+// S_ext_n[LEVEL_N-1]   : Normalized intensity from external source for J -> J'=J-1, S_ext_JJ' / FJJ' = S_ext_n[J']
+//                        S_ext_n = (1 - exp(-TAU_ext)) / (exp(hv/kT_ext) - 1)
 // T                    : Temperature of the cloud (K)
-int coeff_cal(const double *energy_level, double *v, double *E, double *F, double *Br_n, double T) {
+int coeff_cal(const double *energy_level, double *v, double *E, double *F, double *Br_n, double *S_ext_n, double T) {
 	int i;
 	int j, jp;
 	double x;
@@ -121,8 +123,25 @@ int coeff_cal(const double *energy_level, double *v, double *E, double *F, doubl
 	}
 	printf("\n");
 #endif
-	printf("Finished coefficient calculation.\n");
 
+	// S_ext_n[] normalized radiation from external source ========================
+	// S_ext_n = (1 - exp(-TAU_ext)) / (exp(hv/kT_ext) - 1) 
+	// from Planck's law of black-body radiation: B(T,v) = 2hv^3/c^2 / (exp(hv/kT) - 1)
+#if EXT_SOURCE
+#if USE_E_LEVEL_FOR_FREQUENCY
+	x = LIGHT_SPEED * 100.0 * h_CONST / k_CONST / TEMP_EXT;
+	for (j = 0; j < (LEVEL_N - 1); j++) {
+		S_ext_n[j] = (1 - exp(-TAU_EXT)) / (exp(x*(energy_level[j + 1] - energy_level[j])) - 1); //already divided by F = 2hv^3/c^2
+	}
+#else
+	x = h_CONST / k_CONST / TEMP_EXT * 1E9; //the unit of frequency v[] is GHz, GHz = 10^9 Hz
+	for (j = 0; j < (LEVEL_N - 1); j++) {
+		S_ext_n[j] = (1 - exp(-TAU_EXT)) / (exp(x*v[j]) - 1); //already divided by F = 2hv^3/c^2
+	}
+#endif
+#endif
+
+	printf("Finished coefficient calculation.\n");
 	return 0;
 }
 

Codes/coefficients_calculator.h

@@ -4,5 +4,7 @@
 // E[TRANS_N]           : Boltzmann factor, exp(-dEJJ'/kT), for energy difference between J and J', dEJJ'. EJJ' = E[(J-1)J/2+J']
 // F[LEVEL_N-1]         : Flux normalization factor, 2h(vJJ')^3/c^2, FJJ' = F[J']
 // Br_n[LEVEL_N-1]      : Normalized cosmic blackbody radiation intensity for J -> J'=J-1, Br_JJ' = Br[J']
+// S_ext_n[LEVEL_N-1]   : Normalized intensity from external source for J -> J'=J-1, S_ext_JJ' / FJJ' = S_ext_n[J']
+//                        S_ext_n = (1 - exp(-TAU_ext)) / (exp(hv/kT_ext) - 1)
 // T                    : Temperature of the cloud (K)
-int coeff_cal(const double *energy_level, double *v, double *E, double *F, double *Br_n, double T);
+int coeff_cal(const double *energy_level, double *v, double *E, double *F, double *Br_n, double *S_ext_n, double T);

Codes/global_value.h

@@ -6,6 +6,8 @@ extern double E[];            //exp(-hv/kT)
 extern double F[];            //2h(vJJ')^3/c^2, FJJ' = F[J'] (J-1 = J')
 extern double v[];            //frequency (GHz) from J to J'(vJJ'), vJJ' = v[J'] (J-1 = J'), GHz = 10^9 Hz
 extern double Br_n[];         //normalized Cosmic Blackbody Radiation intensity in each level frequency
+extern double S_ext_n[];      // Normalized intensity from external source for J -> J'=J-1, S_ext_JJ' / FJJ' = S_ext_n[J']
+							  // S_ext_n = (1 - exp(-TAU_ext)) / (exp(hv/kT_ext) - 1)
 extern double energy_level[]; //Energy from ground state(J=0) for each level(J)
 extern double a_matrix_i[];   //initial a_matrix[], use 1D array to simulate 2D matrix
 extern double a_matrix[];     //use 1D array to simulate 2D matrix

Codes/main.c

@@ -32,7 +32,9 @@ double C[TRANS_N] = {0.0};            // Collisional excitation rates for J -> J
 double E[TRANS_N] = {0.0};            // Boltzmann factor, exp(-dEJJ'/kT), for energy difference dE between J and J', EJJ' = E[(J-1)J/2+J']
 double F[LEVEL_N-1] = {0.0};          // Flux normalization factor, 2h(vJJ')^3/c^2, FJJ' = F[J']
 double v[LEVEL_N-1] = {0.0};          // Frequency (GHz) for J -> J'=J-1, vJJ' = v[J'], GHz = 10^9 Hz
-double Br_n[LEVEL_N-1] = {0.0};       // Normalized cosmic blackbody radiation intensity for J -> J'=J-1, Br_JJ' = Br[J']
+double Br_n[LEVEL_N-1] = {0.0};       // Normalized cosmic blackbody radiation intensity for J -> J'=J-1, Br_JJ' / FJJ' = Br_n[J']
+double S_ext_n[LEVEL_N-1] = {0.0};    // Normalized intensity from external source for J -> J'=J-1, S_ext_JJ' / FJJ' = S_ext_n[J']
+									  // S_ext_n = (1 - exp(-TAU_ext)) / (exp(hv/kT_ext) - 1)
 double energy_level[LEVEL_N] = {0.0}; // Potential energy (cm^-1) at level J (energy_level[J=0] = 0) (in unit of the inverse of wavelength(l), 1/l)
 double T;                             // Temperature of the cloud (K)
 
@@ -138,7 +140,7 @@ int main()
 #endif	
 
 	// Calculate Coefficients ______________________________________________
-	coeff_cal(energy_level, v, E, F, Br_n, T);		
+	coeff_cal(energy_level, v, E, F, Br_n, S_ext_n, T);
 
 	// Initialize a_matrix_i[] _____________________________________________
 	a_matrix_initialize(a_matrix_i); //fill a_matrix_i[] with C[]
@@ -522,8 +524,10 @@ void testing() {
 
 #if TEST_RATE_EQ_FILL && LEVEL_N == 3
 	test_rate_eq_fill_A_3();       // Pass on 2018.03.12
+#if !EXT_SOURCE
 	test_rate_eq_fill_iso_3();	   // Pass on 2018.03.19
 #endif
+#endif
 
 #if TEST_S_ISO
 	test_source_f_n_iso();         // Pass on 2018.03.14
@@ -541,7 +545,7 @@ void testing() {
 	test_beta_f();                 // Pass on 2018.03.17
 #endif
 
-#if TEST_R_CAL_ISO
+#if TEST_R_CAL_ISO && !EXT_SOURCE
 	test_R_cal_iso();              // Pass on 2018.03.17
 #endif
 

Codes/parameters.h

@@ -21,7 +21,7 @@
 #define WINDOWS                      // Pause at the end of the program and every error
 
 // Number of levels
-#define LEVEL_N 3                       // Total number of levels N (Ex: 2: J = 0 ~ 1)
+#define LEVEL_N 11                       // Total number of levels N (Ex: 2: J = 0 ~ 1)
 #define TOTAL_N ((LEVEL_N+1)*LEVEL_N)/2 // Total number of independent sublevels (size of n[])
 #define TRANS_N ((LEVEL_N-1)*LEVEL_N)/2 // Total transition number for different J (size of C[])
 
@@ -71,12 +71,18 @@
 //#define TAU_INC_RATIO 0.9120108393559098 // Optical depth increase ratio bewteen each step (decreasing)
 
 // Temperature
-#define TEMP_SELE 8                 // The collisional temperature column in the LAMDA data (start from 1) // Before 11
+#define TEMP_SELE 8                  // The collisional temperature column in the LAMDA data (start from 1) // Before 11
 //#define TEMP_SELE 10                 // sio 100K
 //#define TEMP_SELE 31                 // sio 500K
 #define TEMP_B 2.725                 // Cosmic background radiation temperature (K, data from Wiki 2009.12)
 //#define TEMP_B 50.0                  // Use other background radiation temperature (K)
 
+// External continuum point source
+#define EXT_SOURCE 1                 // If there is an external source or not? 0: No, 1: Yes
+#define TEMP_EXT 70.0                // Temperature of the external source (K)
+#define TAU_EXT 1.5                  // Optical depth of the external source
+#define EXT_ANG 0.304*M_PI           // Angle between incident direction of the point source and z-axis
+
 // Integral methods
 #define GSL_INTEGRAL_QNG 0           // QNG non-adaptive Gauss-Kronrod integration
 #define GSL_INTEGRAL_CQUAD 0         // CQUAD doubly-adaptive integration method

Codes/physics_function.c

@@ -189,13 +189,18 @@ void R_cal(const double n[TOTAL_N], double tau, double R[LEVEL_N-1][2]) {
 	// R[j][dM] for the transition from J to j = J' = J-1 with |M - M'| = dM
 	int j; // j = J' = J-1
 	double error;
+	double dRf0, dRf1;
 	unsigned int intervals;
 	Fn_Param params;
 
 	params.n = n;
 	params.tau = tau; //the total optical depth tau[0] in this case(computation)
 	params.k0 = k_f_n(n, cos(TAU_ANG), 0, 0); //[2010.01.22] OBS_ANG was replaced by TAU_ANG
-
+#if EXT_SOURCE
+	dRf0 = 3 * sin(EXT_ANG)*sin(EXT_ANG) / (8 * M_PI);
+	dRf1 = 3 * (1 + cos(EXT_ANG)*cos(EXT_ANG)) / (16 * M_PI);
+#endif
+
 	for (j = 0; j < (LEVEL_N - 1); j++) { //j = 0 ~ LEVEL_N-2
 		params.j = j;
 		R[j][0] = integral(i_f_r0m, &params, &error, &intervals) * 3.0 / 2.0;
@@ -205,7 +210,12 @@ void R_cal(const double n[TOTAL_N], double tau, double R[LEVEL_N-1][2]) {
 
 		R[j][1] = integral(i_f_r1m, &params, &error, &intervals) * 3.0 / 4.0;
 		// R has already integrated over azimuth angle. Thus, it gains a 1/2 factor
-		interval_count += (unsigned long long)intervals;		
+		interval_count += (unsigned long long)intervals;
+
+#if EXT_SOURCE
+		R[j][0] += dRf0 * beta_f(tau_f(n, params.k0, tau, cos(EXT_ANG), 0, j)) * S_ext_n[j];
+		R[j][1] += dRf1 * beta_f(tau_f(n, params.k0, tau, cos(EXT_ANG), 1, j)) * S_ext_n[j];
+#endif
 	}
 }
 
@@ -250,6 +260,24 @@ void output_integrand_R(const double n[TOTAL_N], double tau, int n_sample, const
 	}
 }
 
+// Original way of computing tau function (slow)
+double tau_f(const double n[TOTAL_N], double k0, double tau0, double cos_0, int q, int j) {
+	double cos2 = cos_0 * cos_0;
+	double L_factor, sin2;
+	sin2 = fabs(1 - cos2);  //sin^2(\theta)
+
+#if TwoD //use s*s for 2D velocity field, or c*c for 1D velocity field
+	L_factor = sin(TAU_ANG)*sin(TAU_ANG) / sin2;  // Can be +inf
+#elif OneD
+	L_factor = cos(TAU_ANG)*cos(TAU_ANG) / cos2;  // Can be +inf
+#elif Mix
+	L_factor = (cos(TAU_ANG)*cos(TAU_ANG) + MixRatio * sin(TAU_ANG)*sin(TAU_ANG)) / (cos2 + MixRatio * sin2);
+#else //Isotropic
+	L_factor = 1.;
+#endif
+	return tau0 * (k_f_n(n, cos_0, q, j) / k0) * L_factor;
+}
+
 // Integrand of R_JJ' for dM = 0 (R[J'][0]) without the factor 3/2
 // Return sin^2 * I_pa_n[0]
 double i_f_r0m(double x, void * params) {

Codes/physics_function.h

@@ -46,6 +46,9 @@ void R_cal(const double n[TOTAL_N], double tau, double R[LEVEL_N-1][2]);
 // Given n[], tau[0], number of sample points n_sample, write integrand of R to file output_filename
 void output_integrand_R(const double n[TOTAL_N], double tau, int n_sample, const char* output_filename);
 
+// Original way of computing tau function (slow)
+double tau_f(const double n[TOTAL_N], double k0, double tau0, double cos_0, int q, int j);
+
 // Integrand of R_JJ' for dM = 0 (R[J'][0]) without the factor 3/2
 // Return sin^2 * I_pa_n[0]
 double i_f_r0m(double x, void * params);