From 39e4b71c755bb11cba96a4fe06b67b78708dfe8e Mon Sep 17 00:00:00 2001 From: louis925 Date: Tue, 20 Mar 2018 17:42:26 +0900 Subject: [PATCH] Add option for external background point source --- .gitignore | 5 ++++- Codes/coefficients_calculator.c | 23 +++++++++++++++++++++-- Codes/coefficients_calculator.h | 4 +++- Codes/global_value.h | 2 ++ Codes/main.c | 10 +++++++--- Codes/parameters.h | 10 ++++++++-- Codes/physics_function.c | 32 ++++++++++++++++++++++++++++++-- Codes/physics_function.h | 3 +++ 8 files changed, 78 insertions(+), 11 deletions(-) diff --git a/.gitignore b/.gitignore index 103c4aa..0bc0d83 100644 --- a/.gitignore +++ b/.gitignore @@ -330,4 +330,7 @@ Codes/*.csv Codes/*.dat Codes/*.dem Codes/*.txt -Codes/*.xlsx \ No newline at end of file +Codes/*.xlsx +Codes/*.zip +Codes/*.7z +Codes/*.rar diff --git a/Codes/coefficients_calculator.c b/Codes/coefficients_calculator.c index f2d9411..5affeb6 100644 --- a/Codes/coefficients_calculator.c +++ b/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; } diff --git a/Codes/coefficients_calculator.h b/Codes/coefficients_calculator.h index 7f8d04b..bb3ece3 100644 --- a/Codes/coefficients_calculator.h +++ b/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); diff --git a/Codes/global_value.h b/Codes/global_value.h index a96a448..06d2489 100644 --- a/Codes/global_value.h +++ b/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 diff --git a/Codes/main.c b/Codes/main.c index 1be707b..72c0707 100644 --- a/Codes/main.c +++ b/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 diff --git a/Codes/parameters.h b/Codes/parameters.h index ecd7295..8553624 100644 --- a/Codes/parameters.h +++ b/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 diff --git a/Codes/physics_function.c b/Codes/physics_function.c index 773e9b7..3284247 100644 --- a/Codes/physics_function.c +++ b/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, ¶ms, &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, ¶ms, &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) { diff --git a/Codes/physics_function.h b/Codes/physics_function.h index d3c013e..32993b1 100644 --- a/Codes/physics_function.h +++ b/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);