From 05711da05b3862677403e1af3225ea4fc2d5b33c Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Fri, 21 Mar 2014 13:32:49 +0000 Subject: [PATCH 1/2] Converted frequencies to rest frame of cell in radiation.c --- radiation.c | 32 +++++++++++++++++++++++++++++++- version.h | 2 +- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/radiation.c b/radiation.c index bdc7edbc7..1a1631f31 100755 --- a/radiation.c +++ b/radiation.c @@ -132,12 +132,41 @@ radiation (p, ds) int nconf; double weight_of_packet, y; int ii, jj; + double v_inner[3], v_outer[3], v1, v2; + double freq_inner, freq_outer; + struct photon phot; ii = jj = 0; /* NSH 130605 to remove o3 compile error */ one = &wmain[p->grid]; /* So one is the grid cell of interest */ xplasma = &plasmamain[one->nplasma]; check_plasma (xplasma, "radiation"); - freq = p->freq; + + /* JM 140321 -- #73 Bugfix + In previous version we were comparing these frequencies in the global rest fram + The photon frequency should be shifted to the rest frame of the cell in question + We currently take the average of this frequency along ds. In principle + this could be improved, so we throw an error if the difference between v1 and v2 is large */ + + /* calculate velocity at original position */ + vwind_xyz (p, v_inner); // get velocity vector at new pos + v1 = dot (p->lmn, v_inner); // get direction cosine + + /* Create phot, a photon at the position we are moving to + note that the actual movement of the photon gets done after the call to radiation */ + stuff_phot (p, &phot); // copy photon ptr + move_phot (&phot, ds); // move it by ds + vwind_xyz (&phot, v_outer); // get velocity vector at new pos + v2 = dot (phot.lmn, v_outer); // get direction cosine + + /* calculate photon frequencies in rest frame of cell */ + freq_inner = p->freq * (1. - v1 / C); + freq_outer = phot.freq * (1. - v2 / C); + + /* take the average of the frequencies at original position and original+ds */ + freq = 0.5 * (freq_inner + freq_outer); + + + kappa_tot = frac_ff = kappa_ff (xplasma, freq); /* Add ff opacity */ kappa_tot += frac_comp = kappa_comp (xplasma, freq); /* 70 NSH 1108 calculate compton opacity, store it in kappa_comp and also add it to kappa_tot, the total opacity for the photon path */ @@ -145,6 +174,7 @@ radiation (p, ds) frac_tot = frac_z = 0; /* 59a - ksl - Moved this line out of loop to avoid warning, but notes indicate this is all disagnostic and might be removed */ + if (freq > phot_freq_min) { diff --git a/version.h b/version.h index ab3c3210d..261f331a3 100755 --- a/version.h +++ b/version.h @@ -1,2 +1,2 @@ -#define VERSION "77a_dev" +#define VERSION "ff" #define CHOICE 1 // Compress plasma as much as possible From 9602e55c8c116881e3bc715eaa78a64952bf9960 Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Fri, 21 Mar 2014 16:01:49 +0000 Subject: [PATCH 2/2] Added some comments --- radiation.c | 6 ++++-- version.h | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/radiation.c b/radiation.c index 1a1631f31..64e02b562 100755 --- a/radiation.c +++ b/radiation.c @@ -166,13 +166,15 @@ radiation (p, ds) freq = 0.5 * (freq_inner + freq_outer); - + + /* calculate free-free, compton and ind-compton opacities + note that we also call these with the average frequency along ds */ kappa_tot = frac_ff = kappa_ff (xplasma, freq); /* Add ff opacity */ kappa_tot += frac_comp = kappa_comp (xplasma, freq); /* 70 NSH 1108 calculate compton opacity, store it in kappa_comp and also add it to kappa_tot, the total opacity for the photon path */ kappa_tot += frac_ind_comp = kappa_ind_comp (xplasma, freq); frac_tot = frac_z = 0; /* 59a - ksl - Moved this line out of loop to avoid warning, but notes - indicate this is all disagnostic and might be removed */ + indicate this is all disagnostic and might be removed */ if (freq > phot_freq_min) diff --git a/version.h b/version.h index 261f331a3..ab3c3210d 100755 --- a/version.h +++ b/version.h @@ -1,2 +1,2 @@ -#define VERSION "ff" +#define VERSION "77a_dev" #define CHOICE 1 // Compress plasma as much as possible