From db489b699398dbd0e3ec33dace5853d85934465e Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Thu, 8 Feb 2018 11:14:44 -0600 Subject: [PATCH] Additional diagnostics to track photons as they traverse the grid to address the fu import problem #309. With the -d option is now possible to trance the position of photons at various points in the transfer, most particulary whenever an out of grid warning is generated. --- source/cylindrical.c | 570 ++++++++++++++------------ source/diag.c | 53 ++- source/photon2d.c | 50 ++- source/trans_phot.c | 947 +++++++++++++++++++++++-------------------- 4 files changed, 892 insertions(+), 728 deletions(-) diff --git a/source/cylindrical.c b/source/cylindrical.c index f2560e432..59ee238ce 100755 --- a/source/cylindrical.c +++ b/source/cylindrical.c @@ -40,8 +40,8 @@ double -cylind_ds_in_cell (ndom,p) - int ndom; +cylind_ds_in_cell (ndom, p) + int ndom; PhotPtr p; @@ -62,14 +62,19 @@ cylind_ds_in_cell (ndom,p) if ((p->grid = n = where_in_grid (ndom, p->x)) < 0) - { - Error ("cylind_ds_in_cell: Photon not in grid when routine entered\n"); - return (n); /* Photon was not in wind */ - } + { + if (modes.save_photons) + { + save_photons (p, "NotInGrid_cylind_ds_in_cell"); + } + +//OLD Error ("cylind_ds_in_cell: Photon not in grid when routine entered\n"); + return (n); /* Photon was not in wind */ + } - wind_n_to_ij (ndom, n, &ix, &iz); /*Convert the index n to two dimensions */ + wind_n_to_ij (ndom, n, &ix, &iz); /*Convert the index n to two dimensions */ - smax = VERY_BIG; //initialize smax to a large number + smax = VERY_BIG; //initialize smax to a large number /* Set up the quadratic equations in the radial rho direction */ @@ -77,13 +82,16 @@ cylind_ds_in_cell (ndom,p) b = 2. * (p->lmn[0] * p->x[0] + p->lmn[1] * p->x[1]); c = p->x[0] * p->x[0] + p->x[1] * p->x[1]; - iroot = quadratic (a, b, c - zdom[ndom].wind_x[ix] * zdom[ndom].wind_x[ix], root); /* iroot will be the smallest positive root - if one exists or negative otherwise */ + iroot = quadratic (a, b, c - zdom[ndom].wind_x[ix] * zdom[ndom].wind_x[ix], root); /* iroot will be the smallest positive root + if one exists or negative otherwise */ if (iroot >= 0 && root[iroot] < smax) smax = root[iroot]; - iroot = quadratic (a, b, c - zdom[ndom].wind_x[ix + 1] * zdom[ndom].wind_x[ix + 1], root); + iroot = + quadratic (a, b, + c - zdom[ndom].wind_x[ix + 1] * zdom[ndom].wind_x[ix + 1], + root); if (iroot >= 0 && root[iroot] < smax) smax = root[iroot]; @@ -94,21 +102,21 @@ cylind_ds_in_cell (ndom,p) z1 = zdom[ndom].wind_z[iz]; z2 = zdom[ndom].wind_z[iz + 1]; if (p->x[2] < 0) - { /* We need to worry about which side of the plane the photon is on! */ - z1 *= (-1.); - z2 *= (-1.); - } + { /* We need to worry about which side of the plane the photon is on! */ + z1 *= (-1.); + z2 *= (-1.); + } if (p->lmn[2] != 0.0) - { - q = (z1 - p->x[2]) / p->lmn[2]; - if (q > 0 && q < smax) - smax = q; - q = (z2 - p->x[2]) / p->lmn[2]; - if (q > 0 && q < smax) - smax = q; + { + q = (z1 - p->x[2]) / p->lmn[2]; + if (q > 0 && q < smax) + smax = q; + q = (z2 - p->x[2]) / p->lmn[2]; + if (q > 0 && q < smax) + smax = q; - } + } return (smax); } @@ -153,10 +161,10 @@ cylind_make_grid (ndom, w) one_dom = &zdom[ndom]; if (zdom[ndom].zmax == 0) - { - /* Check if zmax has been set, and if not set it to rmax */ - zdom[ndom].zmax = zdom[ndom].rmax; - } + { + /* Check if zmax has been set, and if not set it to rmax */ + zdom[ndom].zmax = zdom[ndom].rmax; + } Log ("Making cylindrical grid %d\n", ndom); @@ -170,71 +178,85 @@ cylind_make_grid (ndom, w) /* First calculate parameters that are to be calculated at the edge of the grid cell. This is mainly the positions and the velocity */ for (i = 0; i < one_dom->ndim; i++) - { - for (j = 0; j < one_dom->mdim; j++) { - wind_ij_to_n (ndom, i, j, &n); - w[n].x[1] = w[n].xcen[1] = 0; //The cells are all defined in the xz plane - - /*Define the grid points */ - if (one_dom->log_linear == 1) - { // linear intervals - - dr = one_dom->rmax / (one_dom->ndim - 3); - dz = one_dom->zmax / (one_dom->mdim - 3); - w[n].x[0] = i * dr; /* The first zone is at the inner radius of - the wind */ - w[n].x[2] = j * dz; - w[n].xcen[0] = w[n].x[0] + 0.5 * dr; - w[n].xcen[2] = w[n].x[2] + 0.5 * dz; - - } - else - { //logarithmic intervals - - dlogr = (log10 (one_dom->rmax / one_dom->xlog_scale)) / (one_dom->ndim - 3); - dlogz = (log10 (one_dom->zmax / one_dom->zlog_scale)) / (one_dom->mdim - 3); - - if (dlogr <= 0) - { - Error ("cylindrical: dlogr %g is less than 0. This is certainly wrong! Aborting\n", dlogr); - exit (0); - } - - if (dlogz <= 0) - { - Error ("cylindrical: dlogz %g is less than 0. This is certainly wrong! Aborting\n", dlogz); - exit (0); - } - - if (i == 0) - { - w[n].x[0] = 0.0; - w[n].xcen[0] = 0.5 * one_dom->xlog_scale; - } - else - { - w[n].x[0] = one_dom->xlog_scale * pow (10., dlogr * (i - 1)); - w[n].xcen[0] = 0.5 * one_dom->xlog_scale * (pow (10., dlogr * (i - 1)) + pow (10., dlogr * (i))); - } - - if (j == 0) - { - w[n].x[2] = 0.0; - w[n].xcen[2] = 0.5 * one_dom->zlog_scale; - } - else - { - w[n].x[2] = one_dom->zlog_scale * pow (10, dlogz * (j - 1)); - w[n].xcen[2] = 0.5 * one_dom->zlog_scale * (pow (10., dlogz * (j - 1)) + pow (10., dlogz * (j))); - } - } - - xfudge = fmin ((w[n].xcen[0] - w[n].x[0]), (w[n].xcen[2] - w[n].x[2])); - w[n].dfudge = XFUDGE * xfudge; - + for (j = 0; j < one_dom->mdim; j++) + { + wind_ij_to_n (ndom, i, j, &n); + w[n].x[1] = w[n].xcen[1] = 0; //The cells are all defined in the xz plane + + /*Define the grid points */ + if (one_dom->log_linear == 1) + { // linear intervals + + dr = one_dom->rmax / (one_dom->ndim - 3); + dz = one_dom->zmax / (one_dom->mdim - 3); + w[n].x[0] = i * dr; /* The first zone is at the inner radius of + the wind */ + w[n].x[2] = j * dz; + w[n].xcen[0] = w[n].x[0] + 0.5 * dr; + w[n].xcen[2] = w[n].x[2] + 0.5 * dz; + + } + else + { //logarithmic intervals + + dlogr = + (log10 (one_dom->rmax / one_dom->xlog_scale)) / + (one_dom->ndim - 3); + dlogz = + (log10 (one_dom->zmax / one_dom->zlog_scale)) / + (one_dom->mdim - 3); + + if (dlogr <= 0) + { + Error + ("cylindrical: dlogr %g is less than 0. This is certainly wrong! Aborting\n", + dlogr); + exit (0); + } + + if (dlogz <= 0) + { + Error + ("cylindrical: dlogz %g is less than 0. This is certainly wrong! Aborting\n", + dlogz); + exit (0); + } + + if (i == 0) + { + w[n].x[0] = 0.0; + w[n].xcen[0] = 0.5 * one_dom->xlog_scale; + } + else + { + w[n].x[0] = + one_dom->xlog_scale * pow (10., dlogr * (i - 1)); + w[n].xcen[0] = + 0.5 * one_dom->xlog_scale * (pow (10., dlogr * (i - 1)) + + pow (10., dlogr * (i))); + } + + if (j == 0) + { + w[n].x[2] = 0.0; + w[n].xcen[2] = 0.5 * one_dom->zlog_scale; + } + else + { + w[n].x[2] = one_dom->zlog_scale * pow (10, dlogz * (j - 1)); + w[n].xcen[2] = + 0.5 * one_dom->zlog_scale * (pow (10., dlogz * (j - 1)) + + pow (10., dlogz * (j))); + } + } + + xfudge = + fmin ((w[n].xcen[0] - w[n].x[0]), (w[n].xcen[2] - w[n].x[2])); + w[n].dfudge = XFUDGE * xfudge; + + } } - } return (0); } @@ -274,14 +296,20 @@ cylind_wind_complete (ndom, w) one_dom->wind_z[j] = w[nstart + j].x[2]; for (i = 0; i < one_dom->ndim - 1; i++) - one_dom->wind_midx[i] = 0.5 * (w[nstart + i * mdim].x[0] + w[nstart + (i + 1) * mdim].x[0]); + one_dom->wind_midx[i] = + 0.5 * (w[nstart + i * mdim].x[0] + w[nstart + (i + 1) * mdim].x[0]); for (j = 0; j < one_dom->mdim - 1; j++) - one_dom->wind_midz[j] = 0.5 * (w[nstart + j].x[2] + w[nstart + j + 1].x[2]); + one_dom->wind_midz[j] = + 0.5 * (w[nstart + j].x[2] + w[nstart + j + 1].x[2]); /* Add something plausible for the edges */ - one_dom->wind_midx[one_dom->ndim - 1] = 2. * one_dom->wind_x[nstart + ndim - 1] - one_dom->wind_midx[nstart + ndim - 2]; - one_dom->wind_midz[one_dom->mdim - 1] = 2. * one_dom->wind_z[nstart + mdim - 1] - one_dom->wind_midz[nstart + mdim - 2]; + one_dom->wind_midx[one_dom->ndim - 1] = + 2. * one_dom->wind_x[nstart + ndim - 1] - one_dom->wind_midx[nstart + + ndim - 2]; + one_dom->wind_midz[one_dom->mdim - 1] = + 2. * one_dom->wind_z[nstart + mdim - 1] - one_dom->wind_midz[nstart + + mdim - 2]; return (0); } @@ -336,7 +364,7 @@ cylind_wind_complete (ndom, w) int cylind_volumes (ndom, w) - int ndom; // domain number + int ndom; // domain number WindPtr w; { int i, j, n; @@ -354,101 +382,106 @@ cylind_volumes (ndom, w) for (i = 0; i < one_dom->ndim; i++) - { - for (j = 0; j < one_dom->mdim; j++) { - - wind_ij_to_n (ndom, i, j, &n); - - rmin = one_dom->wind_x[i]; - rmax = one_dom->wind_x[i + 1]; - zmin = one_dom->wind_z[j]; - zmax = one_dom->wind_z[j + 1]; - - /* this is the full cell volume, which is adjusted by the fraction - of the cell that is in the wind later if necessary */ - //leading factor of 2 added to allow for volume above and below plane (SSMay04) - cell_volume = 2 * PI * (rmax * rmax - rmin * rmin) * (zmax - zmin); - - // XXX Why is it necessary to do the check indicated by the if statement. - /* JM 1711 -- only try to assign the cell if it has not already been assigned */ - if (w[n].inwind == W_NOT_ASSIGNED) - { - if (one_dom->wind_type == IMPORT) { - Error("Shouldn't be redefining inwind in cylind_volumes with imported model.\n"); - exit(0); - } - - n_inwind = cylind_is_cell_in_wind (n); - - if (n_inwind == W_NOT_INWIND) - { - fraction = 0.0; - jj = 0; - kk = RESOLUTION * RESOLUTION; - } - else if (n_inwind == W_ALL_INWIND) - { - fraction = 1.0; - jj = kk = RESOLUTION * RESOLUTION; - } - else - { /* Determine whether the grid cell is in the wind */ - num = denom = 0; - jj = kk = 0; - dr = (rmax - rmin) / RESOLUTION; - dz = (zmax - zmin) / RESOLUTION; - for (r = rmin + dr / 2; r < rmax; r += dr) - { - for (z = zmin + dz / 2; z < zmax; z += dz) - { - denom += r * r; - kk++; - x[0] = r; - x[1] = 0; - x[2] = z; - - if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) - { - num += r * r; /* 0 implies in wind */ - jj++; - } - } - } - fraction = num / denom; - } - - /* OK now make the final assignement of nwind and fix the volumes */ - /* XXX JM - not clear why these additional if statements are necessary */ - if (jj == 0) - { - w[n].inwind = W_NOT_INWIND; // The cell is not in the wind - w[n].vol = 0.0; - } - else if (jj == kk) - { - w[n].inwind = W_ALL_INWIND; // All of cell is inwind - w[n].vol = cell_volume; - } - else - { - w[n].inwind = W_PART_INWIND; // Some of cell is inwind - w[n].vol = cell_volume * fraction; - } - } - - /* JM 1711 -- the following two if statements are for if the inwind values are - already assigned, for example by an imported model */ - /* need to zero volumes for cells not in the wind */ - else if (w[n].inwind == W_NOT_INWIND) { - w[n].vol = 0.0; - } - - else if (w[n].inwind == W_ALL_INWIND) { - w[n].vol = cell_volume; - } + for (j = 0; j < one_dom->mdim; j++) + { + + wind_ij_to_n (ndom, i, j, &n); + + rmin = one_dom->wind_x[i]; + rmax = one_dom->wind_x[i + 1]; + zmin = one_dom->wind_z[j]; + zmax = one_dom->wind_z[j + 1]; + + /* this is the full cell volume, which is adjusted by the fraction + of the cell that is in the wind later if necessary */ + //leading factor of 2 added to allow for volume above and below plane (SSMay04) + cell_volume = 2 * PI * (rmax * rmax - rmin * rmin) * (zmax - zmin); + + // XXX Why is it necessary to do the check indicated by the if statement. + /* JM 1711 -- only try to assign the cell if it has not already been assigned */ + if (w[n].inwind == W_NOT_ASSIGNED) + { + if (one_dom->wind_type == IMPORT) + { + Error + ("Shouldn't be redefining inwind in cylind_volumes with imported model.\n"); + exit (0); + } + + n_inwind = cylind_is_cell_in_wind (n); + + if (n_inwind == W_NOT_INWIND) + { + fraction = 0.0; + jj = 0; + kk = RESOLUTION * RESOLUTION; + } + else if (n_inwind == W_ALL_INWIND) + { + fraction = 1.0; + jj = kk = RESOLUTION * RESOLUTION; + } + else + { /* Determine whether the grid cell is in the wind */ + num = denom = 0; + jj = kk = 0; + dr = (rmax - rmin) / RESOLUTION; + dz = (zmax - zmin) / RESOLUTION; + for (r = rmin + dr / 2; r < rmax; r += dr) + { + for (z = zmin + dz / 2; z < zmax; z += dz) + { + denom += r * r; + kk++; + x[0] = r; + x[1] = 0; + x[2] = z; + + if (where_in_wind (x, &ndomain) == W_ALL_INWIND + && ndom == ndomain) + { + num += r * r; /* 0 implies in wind */ + jj++; + } + } + } + fraction = num / denom; + } + + /* OK now make the final assignement of nwind and fix the volumes */ + /* XXX JM - not clear why these additional if statements are necessary */ + if (jj == 0) + { + w[n].inwind = W_NOT_INWIND; // The cell is not in the wind + w[n].vol = 0.0; + } + else if (jj == kk) + { + w[n].inwind = W_ALL_INWIND; // All of cell is inwind + w[n].vol = cell_volume; + } + else + { + w[n].inwind = W_PART_INWIND; // Some of cell is inwind + w[n].vol = cell_volume * fraction; + } + } + + /* JM 1711 -- the following two if statements are for if the inwind values are + already assigned, for example by an imported model */ + /* need to zero volumes for cells not in the wind */ + else if (w[n].inwind == W_NOT_INWIND) + { + w[n].vol = 0.0; + } + + else if (w[n].inwind == W_ALL_INWIND) + { + w[n].vol = cell_volume; + } + } } - } return (0); } @@ -502,17 +535,18 @@ cylind_where_in_grid (ndom, x) one_dom = &zdom[ndom]; - z = fabs (x[2]); /* This is necessary to get correct answer above - and below plane */ + z = fabs (x[2]); /* This is necessary to get correct answer above + and below plane */ if (z == 0) - z = 1.e4; //Force z to be positive 02feb ksl - rho = sqrt (x[0] * x[0] + x[1] * x[1]); /* This is distance from z - axis */ + z = 1.e4; //Force z to be positive 02feb ksl + rho = sqrt (x[0] * x[0] + x[1] * x[1]); /* This is distance from z + axis */ /* Check to see if x is outside the region of the calculation */ - if (rho > one_dom->wind_x[one_dom->ndim - 1] || z > one_dom->wind_z[one_dom->mdim - 1]) - { - return (-2); /* x is outside grid */ - } + if (rho > one_dom->wind_x[one_dom->ndim - 1] + || z > one_dom->wind_z[one_dom->mdim - 1]) + { + return (-2); /* x is outside grid */ + } if (rho < one_dom->wind_x[0]) return (-1); @@ -558,8 +592,8 @@ cylind_where_in_grid (ndom, x) int cylind_get_random_location (n, x) - int n; // Cell in which to create postion - double x[]; // Returned position + int n; // Cell in which to create postion + double x[]; // Returned position { int i, j; int inwind; @@ -583,25 +617,27 @@ cylind_get_random_location (n, x) /* Generate a position which is both in the cell and in the wind */ inwind = W_NOT_INWIND; while (inwind != W_ALL_INWIND || ndomain != ndom) - { - r = sqrt (rmin * rmin + (rand () / (MAXRAND - 0.5)) * (rmax * rmax - rmin * rmin)); + { + r = + sqrt (rmin * rmin + + (rand () / (MAXRAND - 0.5)) * (rmax * rmax - rmin * rmin)); // Generate the azimuthal location - phi = 2. * PI * (rand () / MAXRAND); - x[0] = r * cos (phi); - x[1] = r * sin (phi); + phi = 2. * PI * (rand () / MAXRAND); + x[0] = r * cos (phi); + x[1] = r * sin (phi); - x[2] = zmin + (zmax - zmin) * (rand () / (MAXRAND - 0.5)); - inwind = where_in_wind (x, &ndomain); /* Some photons will not be in the wind - because the boundaries of the wind split the grid cell */ - } + x[2] = zmin + (zmax - zmin) * (rand () / (MAXRAND - 0.5)); + inwind = where_in_wind (x, &ndomain); /* Some photons will not be in the wind + because the boundaries of the wind split the grid cell */ + } - zz = rand () / MAXRAND - 0.5; //positions above are all at +z distances + zz = rand () / MAXRAND - 0.5; //positions above are all at +z distances if (zz < 0) - x[2] *= -1; /* The photon is in the bottom half of the wind */ + x[2] *= -1; /* The photon is in the bottom half of the wind */ return (inwind); } @@ -661,34 +697,34 @@ cylind_extend_density (ndom, w) */ for (i = 0; i < ndim - 1; i++) - { - for (j = 0; j < mdim - 1; j++) { - wind_ij_to_n (ndom, i, j, &n); - if (w[n].vol == 0) - - { //Then this grid point is not in the wind - - wind_ij_to_n (ndom, i + 1, j, &m); - if (w[m].vol > 0) - { //Then the windcell in the +x direction is in the wind and - // we can copy the densities to the grid cell n - w[n].nplasma = w[m].nplasma; - - } - else if (i > 0) - { - wind_ij_to_n (ndom, i - 1, j, &m); - if (w[m].vol > 0) - { //Then the grid cell in the -x direction is in the wind and - // we can copy the densities to the grid cell n - w[n].nplasma = w[m].nplasma; - - } - } - } + for (j = 0; j < mdim - 1; j++) + { + wind_ij_to_n (ndom, i, j, &n); + if (w[n].vol == 0) + + { //Then this grid point is not in the wind + + wind_ij_to_n (ndom, i + 1, j, &m); + if (w[m].vol > 0) + { //Then the windcell in the +x direction is in the wind and + // we can copy the densities to the grid cell n + w[n].nplasma = w[m].nplasma; + + } + else if (i > 0) + { + wind_ij_to_n (ndom, i - 1, j, &m); + if (w[m].vol > 0) + { //Then the grid cell in the -x direction is in the wind and + // we can copy the densities to the grid cell n + w[n].nplasma = w[m].nplasma; + + } + } + } + } } - } return (0); @@ -718,7 +754,7 @@ Note that it simply calls where_in_wind multiple times. int cylind_is_cell_in_wind (n) - int n; // cell number + int n; // cell number { int i, j; double r, z, dr, dz; @@ -733,18 +769,18 @@ cylind_is_cell_in_wind (n) /* First check if the cell is in the boundary */ if (i >= (one_dom->ndim - 2) && j >= (one_dom->mdim - 2)) - { - return (W_NOT_INWIND); - } + { + return (W_NOT_INWIND); + } /* Assume that if all four corners are in the wind that the entire cell is in the wind. check_corners also now checks that the corners are in the wind of a specific domain */ if (check_corners_inwind (n) == 4) - { - return (W_ALL_INWIND); - } + { + return (W_ALL_INWIND); + } /* So at this point, we have dealt with the easy cases of being in the region of the wind which is used for the @@ -765,46 +801,46 @@ cylind_is_cell_in_wind (n) x[1] = 0; for (z = zmin + dz / 2; z < zmax; z += dz) - { - x[2] = z; + { + x[2] = z; - x[0] = rmin; + x[0] = rmin; - if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) - { - return (W_PART_INWIND); - } + if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) + { + return (W_PART_INWIND); + } - x[0] = rmax; + x[0] = rmax; - if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) - { - return (W_PART_INWIND); + if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) + { + return (W_PART_INWIND); + } } - } // Check inner and outer boundary in the z direction for (r = rmin + dr / 2; r < rmax; r += dr) - { + { - x[0] = r; + x[0] = r; - x[2] = zmin; + x[2] = zmin; - if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) - { - return (W_PART_INWIND); - } + if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) + { + return (W_PART_INWIND); + } - x[2] = zmax; + x[2] = zmax; - if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) - { - return (W_PART_INWIND); + if (where_in_wind (x, &ndomain) == W_ALL_INWIND && ndom == ndomain) + { + return (W_PART_INWIND); + } } - } /* If one has reached this point, then this wind cell is not in the wind */ return (W_NOT_INWIND); diff --git a/source/diag.c b/source/diag.c index a3f66050c..49ec38e3d 100755 --- a/source/diag.c +++ b/source/diag.c @@ -173,7 +173,7 @@ init_extra_diagnostics () if (eplinit == 0 && modes.extra_diagnostics) { - epltptr = fopen ("python.ext", "w"); + epltptr = fopen ("python.ext.txt", "w"); eplinit = 1; } @@ -274,6 +274,28 @@ save_photon_stats (one, p, ds, w_ave) return (0); } + +/*********************************************************** + Space Telescope Science Insittue + +Synopsis: + save_extract_photons saves informations about phtoons in + a particulare wavelength gange + +Arguments: +Returns: + +Description: + +Notes: + Moved here to save duplicating code between bf_estimators_increment and radiation. + +History: + 1802 ksl Functionality moved from extract.c to consolidiate how extra + diagnostics were carried out. + +**************************************************************/ + int save_extract_photons (n, p, pp, v) int n; @@ -289,15 +311,40 @@ save_extract_photons (n, p, pp, v) return (0); } +/*********************************************************** + Space Telescope Science Institute + +Synopsis: + save_photon + +Arguments: + p Photon pointer + comment An arbitrary comment + +Returns: + +Description: + This is a diagnositc that allows one to print out information about + a photon at any time. It was written as part of the effort to + debub imports for the fu_ori project. + +Notes: + +History: + !802 ksl Coded + +**************************************************************/ + + int save_photons (p, comment) PhotPtr p; char comment[]; { fprintf (epltptr, -"PHOTON %3d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %s \n", +"PHOTON %3d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %3d %3d %s \n", p->np, p->x[0], p->x[1], p->x[2], p->lmn[0], p->lmn[1], -p->lmn[2], comment); +p->lmn[2], p->grid, p->istat,comment); return(0); } diff --git a/source/photon2d.c b/source/photon2d.c index d9aca57ae..9070c6cc7 100755 --- a/source/photon2d.c +++ b/source/photon2d.c @@ -58,18 +58,18 @@ translate (w, pp, tau_scat, tau, nres) if (where_in_wind (pp->x, &ndomain) < 0) { istat = translate_in_space (pp); - if (modes.save_photons) - { - save_photons(pp,"Space"); - } + if (modes.save_photons) + { + save_photons (pp, "Space"); + } } else if ((pp->grid = where_in_grid (ndomain, pp->x)) >= 0) { istat = translate_in_wind (w, pp, tau_scat, tau, nres); if (modes.save_photons) - { - save_photons(pp,"Wind"); - } + { + save_photons (pp, "Wind"); + } } else { @@ -128,14 +128,17 @@ translate_in_space (pp) { stuff_phot (pp, &ptest); move_phot (&ptest, ds + DFUDGE); /* So now ptest is at the edge of the wind as defined by the boundary - From here on we should be in the grid */ + From here on we should be in the grid */ /* XXX this is a test. We check at the start whether we are in the grid */ if ((ifail = where_in_grid (ndom, ptest.x)) < 0) { - Error ("translate_in_space: Failure %10.3e %10.3e %10.3e %d %d %d\n", ptest.x[0], ptest.x[1], ptest.x[2],ptest.np,ifail, xxxbound); - }; + if (modes.save_photons) + { + save_photons (pp, "NotInGrid_translate_in_space1"); + } + } /* XXX this ends the test */ @@ -161,10 +164,11 @@ translate_in_space (pp) } else { - Error - ("translate_in_space: Photon not in grid: %10.3e %10.3e %10.3e %03d %10.3e %10.3e\n", - ptest.x[0], ptest.x[1], ptest.x[2],ptest.np,s,smax); - break; + if (modes.save_photons) + { + save_photons (pp, "NotInGrid_translate_in_space2"); + } + break; } } @@ -527,9 +531,10 @@ return and record an error */ { if (translate_in_wind_failure < 1000) { - Error - ("translate_in_wind: Photon not in grid when routine entered\n"); - translate_in_wind_failure += 1; + if (modes.save_photons) + { + save_photons (p, "NotInGrid_translate_in_wind"); + } } return (n); /* Photon was not in grid */ } @@ -578,8 +583,8 @@ return and record an error */ Error ("translate_in_wind: Grid cell %d of photon is not in wind, moving photon %.2e\n", n, smax); - Error ("translate_in_wind: photon %d position: x %g y %g z %g\n", p->np, - p->x[0], p->x[1], p->x[2]); + Error ("translate_in_wind: photon %d position: x %g y %g z %g\n", + p->np, p->x[0], p->x[1], p->x[2]); move_phot (p, smax); return (p->istat); @@ -707,8 +712,11 @@ return and record an error */ if ((p->grid = n = where_in_grid (ndom, p->x)) < 0) { - Error ("ds_in_cell: Photon not in grid when routine entered\n"); - return (n); /* Photon was not in grid */ + if (modes.save_photons) + { + save_photons (p, "NotInGrid_ds_in_cell"); + } + return(n); } /* Assign the pointers for the cell containing the photon */ diff --git a/source/trans_phot.c b/source/trans_phot.c index 9f4b9ca85..56342424b 100755 --- a/source/trans_phot.c +++ b/source/trans_phot.c @@ -16,6 +16,9 @@ !0 -> the normal option for python and hence the need to call "extract" Returns: + + At the end of the routine p is the last point where the photon was in the wind, + not the outer boundary of the radiative transfer Description: This routine oversees the propagation of individual photons. The main loop @@ -58,13 +61,15 @@ FILE *pltptr; int plinit = 0; long n_lost_to_dfudge = 0; -int trans_phot (WindPtr w, PhotPtr p, int iextract /* 0 means do not extract along specific angles; nonzero implies to extract */ - ) +/* 0 means do not extract along specific angles; nonzero implies to extract */ + +int +trans_phot (WindPtr w, PhotPtr p, int iextract) { int nphot; struct photon pp, pextract; int nnscat; - int absorb_reflect; /* this is a variable used to store geo.absorb_reflect during exxtract */ + int absorb_reflect; /* this is a variable used to store geo.absorb_reflect during exxtract */ int nerr; double p_norm, tau_norm; @@ -72,147 +77,171 @@ int trans_phot (WindPtr w, PhotPtr p, int iextract /* 0 means do not extrac /* 05jul -- not clear whether this is needed and why it is different from DEBUG */ /* 1411 -- JM -- Debug usage has been altered. See #111, #120 */ + /* 1802 XXX ksl - This needs to be incoropratred into diag.c */ if (modes.track_resonant_scatters && plinit == 0) - { - pltptr = fopen ("python.xyz", "w"); - plinit = 1; - } + { + pltptr = fopen ("python.xyz", "w"); + plinit = 1; + } Log ("\n"); for (nphot = 0; nphot < NPHOT; nphot++) - { - - // This is just a watchdog method to tell the user the program is still running - - if (nphot % 50000 == 0) - Log ("Photon %7d of %7d or %6.3f per cent \n", nphot, NPHOT, nphot * 100. / NPHOT); - - Log_flush (); - - /* 74a_ksl Check that the weights are real */ - - if (sane_check (p[nphot].w)) { - Error ("trans_phot:sane_check photon %d has weight %e\n", nphot, p[nphot].w); - } - /* Next block added by SS Jan 05 - for anisotropic scattering with extract we want to be sure that everything is - initialised (by scatter?) before calling extract for macro atom photons. Insert this call to scatter which should do - this. */ + // This is just a watchdog method to tell the user the program is still running + + if (nphot % 50000 == 0) + Log ("Photon %7d of %7d or %6.3f per cent \n", nphot, NPHOT, + nphot * 100. / NPHOT); + + Log_flush (); + + /* 74a_ksl Check that the weights are real */ + + if (sane_check (p[nphot].w)) + { + Error ("trans_phot:sane_check photon %d has weight %e\n", nphot, + p[nphot].w); + } + /* Next block added by SS Jan 05 - for anisotropic scattering with extract we want to be sure that everything is + initialised (by scatter?) before calling extract for macro atom photons. Insert this call to scatter which should do + this. */ + + + if (geo.rt_mode == RT_MODE_MACRO + && geo.scatter_mode == SCATTER_MODE_ANISOTROPIC) + { + if (p[nphot].origin == PTYPE_WIND) + { + if (p[nphot].nres > -1 && p[nphot].nres < NLINES) + { + geo.rt_mode = RT_MODE_2LEVEL; + /* 74a_ksl Check to see when a photon weight is becoming unreal */ + if (sane_check (p[nphot].w)) + { + Error + ("trans_phot:sane_check photon %d has weight %e before scatter\n", + nphot, p[nphot].w); + } + if ((nerr = + scatter (&p[nphot], &p[nphot].nres, &nnscat)) != 0) + { + Error + ("trans_phot: Bad return from scatter %d at point 1", + nerr); + } + /* 74a_ksl Check to see when a photon weight is becoming unreal */ + if (sane_check (p[nphot].w)) + { + Error + ("trans_phot:sane_check photon %d has weight %e aftger scatter\n", + nphot, p[nphot].w); + } + geo.rt_mode = RT_MODE_MACRO; + } + } + } + + + stuff_phot (&p[nphot], &pp); + absorb_reflect = geo.absorb_reflect; + + /* The next if statement is executed if we are calculating the detailed spectrum and makes sure we always run extract on + the original photon no matter where it was generated */ - if (geo.rt_mode == RT_MODE_MACRO && geo.scatter_mode == SCATTER_MODE_ANISOTROPIC) - { - if (p[nphot].origin == PTYPE_WIND) - { - if (p[nphot].nres > -1 && p[nphot].nres < NLINES) - { - geo.rt_mode = RT_MODE_2LEVEL; - /* 74a_ksl Check to see when a photon weight is becoming unreal */ - if (sane_check (p[nphot].w)) - { - Error ("trans_phot:sane_check photon %d has weight %e before scatter\n", nphot, p[nphot].w); - } - if ((nerr = scatter (&p[nphot], &p[nphot].nres, &nnscat)) != 0) - { - Error ("trans_phot: Bad return from scatter %d at point 1", nerr); - } - /* 74a_ksl Check to see when a photon weight is becoming unreal */ - if (sane_check (p[nphot].w)) - { - Error ("trans_phot:sane_check photon %d has weight %e aftger scatter\n", nphot, p[nphot].w); - } - geo.rt_mode = RT_MODE_MACRO; - } - } - } + if (iextract) + { + // SS - for reflecting disk have to make sure disk photons are only extracted once. Note we restore the + // correct geo.absorb_reflect value as soon as the photons are extracted! + if (absorb_reflect == BACK_RAD_SCATTER + && p[nphot].origin == PTYPE_DISK) + { + geo.absorb_reflect = BACK_RAD_ABSORB_AND_DESTROY; + } - stuff_phot (&p[nphot], &pp); - absorb_reflect = geo.absorb_reflect; - /* The next if statement is executed if we are calculating the detailed spectrum and makes sure we always run extract on - the original photon no matter where it was generated */ + stuff_phot (&p[nphot], &pextract); - if (iextract) - { - // SS - for reflecting disk have to make sure disk photons are only extracted once. Note we restore the - // correct geo.absorb_reflect value as soon as the photons are extracted! - if (absorb_reflect == BACK_RAD_SCATTER && p[nphot].origin == PTYPE_DISK) - { - geo.absorb_reflect = BACK_RAD_ABSORB_AND_DESTROY; - } + /* We then increase weight to account for number of scatters. This is done because in extract we multiply by the escape + probability along a given direction, but we also need to divide the weight by the mean escape probability, which is + equal to 1/nnscat */ + if (geo.scatter_mode == SCATTER_MODE_THERMAL + && pextract.nres <= NLINES && pextract.nres > 0) + { + /* we normalised our rejection method by the escape probability along the vector of maximum velocity gradient. + First find the sobolev optical depth along that vector */ + tau_norm = + sobolev (&wmain[pextract.grid], pextract.x, -1.0, + lin_ptr[pextract.nres], + wmain[pextract.grid].dvds_max); - stuff_phot (&p[nphot], &pextract); + /* then turn into a probability */ + p_norm = p_escape_from_tau (tau_norm); + } + else + { + p_norm = 1.0; - /* We then increase weight to account for number of scatters. This is done because in extract we multiply by the escape - probability along a given direction, but we also need to divide the weight by the mean escape probability, which is - equal to 1/nnscat */ - if (geo.scatter_mode == SCATTER_MODE_THERMAL && pextract.nres <= NLINES && pextract.nres > 0) - { - /* we normalised our rejection method by the escape probability along the vector of maximum velocity gradient. - First find the sobolev optical depth along that vector */ + /* throw an error if nnscat does not equal 1 */ + if (pextract.nnscat != 1) + Error + ("trans_phot: nnscat is %i for photon %i in scatter mode %i! nres %i NLINES %i\n", + pextract.nnscat, nphot, geo.scatter_mode, pextract.nres, + NLINES); + } - tau_norm = sobolev (&wmain[pextract.grid], pextract.x, -1.0, lin_ptr[pextract.nres], wmain[pextract.grid].dvds_max); - /* then turn into a probability */ - p_norm = p_escape_from_tau (tau_norm); - } - else - { - p_norm = 1.0; + /* We then increase weight to account for number of scatters. This is done because in extract we multiply by the escape + probability along a given direction, but we also need to divide the weight by the mean escape probability, which is + equal to 1/nnscat */ - /* throw an error if nnscat does not equal 1 */ - if (pextract.nnscat != 1) - Error - ("trans_phot: nnscat is %i for photon %i in scatter mode %i! nres %i NLINES %i\n", - pextract.nnscat, nphot, geo.scatter_mode, pextract.nres, NLINES); - } + pextract.w *= p[nphot].nnscat / p_norm; + if (sane_check (pextract.w)) + { + Error + ("trans_phot: sane_check photon %d has weight %e before extract\n", + nphot, pextract.w); + } + extract (w, &pextract, pextract.origin); - /* We then increase weight to account for number of scatters. This is done because in extract we multiply by the escape - probability along a given direction, but we also need to divide the weight by the mean escape probability, which is - equal to 1/nnscat */ + // Restore the correct disk illumination + if (absorb_reflect == BACK_RAD_SCATTER + && p[nphot].origin == PTYPE_DISK) + { + geo.absorb_reflect = BACK_RAD_SCATTER; + } + } /* End of extract loop */ - pextract.w *= p[nphot].nnscat / p_norm; + p[nphot].np = nphot; + trans_phot_single (w, &p[nphot], iextract); - if (sane_check (pextract.w)) - { - Error ("trans_phot: sane_check photon %d has weight %e before extract\n", nphot, pextract.w); - } - extract (w, &pextract, pextract.origin); - // Restore the correct disk illumination - if (absorb_reflect == BACK_RAD_SCATTER && p[nphot].origin == PTYPE_DISK) - { - geo.absorb_reflect = BACK_RAD_SCATTER; - } } - p[nphot].np = nphot; - trans_phot_single (w, &p[nphot], iextract); - - } - /* This is the end of the loop over all of the photons; after this the routine returns */ - + // 130624 ksl Line added to complete watchdog timer, Log ("\n\n"); /* sometimes photons scatter near the edge of the wind and get pushed out by DFUDGE. We record these */ if (n_lost_to_dfudge > 0) - Error ("trans_phot: %ld photons were lost due to DFUDGE (=%8.4e) pushing them outside of the wind after scatter\n", n_lost_to_dfudge, DFUDGE); + Error + ("trans_phot: %ld photons were lost due to DFUDGE (=%8.4e) pushing them outside of the wind after scatter\n", + n_lost_to_dfudge, DFUDGE); - n_lost_to_dfudge = 0; // reset the counter + n_lost_to_dfudge = 0; // reset the counter return (0); } @@ -273,351 +302,395 @@ trans_phot_single (WindPtr w, PhotPtr p, int iextract) tau = 0; icell = 0; - n = 0; // Needed to avoid 03 warning, but it is not clear that it is defined as expected. + n = 0; // Needed to avoid 03 warning, but it is not clear that it is defined as expected. - /* This is the beginning of the loop for each photon and executes until the photon leaves the wind */ + /* This is the beginning of the loop for a single photon and executes until the photon leaves the wind */ while (istat == P_INWIND) - { - - /* translate involves only a single shell (or alternatively a single tranfer in the windless region). istat as returned by - should either be 0 in which case the photon hit the other side of the shell without scattering or 1 in which case there - was a scattering event in the shell, 2 in which case the photon reached the outside edge of the grid and escaped, 3 in - which case it reach the inner edge and was reabsorbed. If the photon escapes then we leave the photon at the position - of it's last scatter. In most other cases though we store the final position of the photon. */ - - - istat = translate (w, &pp, tau_scat, &tau, &nres); - /* nres is the resonance at which the photon was stopped. At present the same value is also stored in pp->nres, but I have - not yet eliminated it from translate. ?? 02jan ksl */ - - icell++; - istat = walls (&pp, p,normal); - // pp is where the photon is going, p is where it was - - if (modes.ispy) - ispy (&pp, n); - - if (istat == -1) { - Error_silent ("trans_phot: Abnormal return from translate on photon %d\n", p->np); - break; - } - if (pp.w < weight_min) - { - istat = pp.istat = P_ABSORB; /* This photon was absorbed by continuum opacity within the wind */ - pp.tau = VERY_BIG; - stuff_phot (&pp, p); - break; - } + /* translate involves only a single shell (or alternatively a single tranfer in the windless region). istat as returned by + should either be 0 in which case the photon hit the other side of the shell without scattering or 1 in which case there + was a scattering event in the shell, 2 in which case the photon reached the outside edge of the grid and escaped, 3 in + which case it reach the inner edge and was reabsorbed. If the photon escapes then we leave the photon at the position + of it's last scatter. In most other cases though we store the final position of the photon. */ + + + istat = translate (w, &pp, tau_scat, &tau, &nres); + /* nres is the resonance at which the photon was stopped. At present the same value is also stored in pp->nres, but I have + not yet eliminated it from translate. ?? 02jan ksl */ + + icell++; + istat = walls (&pp, p, normal); + // pp is where the photon is going, p is where it was + + if (modes.ispy) + ispy (&pp, n); + + if (istat == -1) + { + Error_silent + ("trans_phot: Abnormal return from translate on photon %d\n", + p->np); + break; + } + + if (pp.w < weight_min) + { + istat = pp.istat = P_ABSORB; /* This photon was absorbed by continuum opacity within the wind */ + pp.tau = VERY_BIG; + stuff_phot (&pp, p); + break; + } + + if (istat == P_HIT_STAR) + { /* It hit the star */ + geo.lum_star_back += pp.w; + if (geo.absorb_reflect == BACK_RAD_SCATTER) + { + randvcos (pp.lmn, normal); + stuff_phot (&pp, p); + tau_scat = -log (1. - (rand () + 0.5) / MAXRAND); + istat = pp.istat = P_INWIND; // if we got here, the photon stays in the wind- make sure istat doesn't say scattered still! + tau = 0; + if (iextract) + { + stuff_phot (&pp, &pextract); + extract (w, &pextract, PTYPE_STAR); // Treat as wind photon for purpose of extraction + } + } + else + { /*This is the end of the line for this photon */ + stuff_phot (&pp, p); + break; + } + } + + if (istat == P_HIT_DISK) + { + /* It hit the disk */ + + /* Store the energy of the photon bundle into a disk structure so that one can determine later how much and where the + disk was heated by photons. + Note that the disk is defined from 0 to NRINGS-2. NRINGS-1 contains the position of the outer radius of the disk. */ + + rrr = sqrt (dot (pp.x, pp.x)); + kkk = 0; + while (rrr > qdisk.r[kkk] && kkk < NRINGS - 1) + kkk++; + kkk--; // So that the heating refers to the heating between kkk and kkk+1 + qdisk.nhit[kkk]++; + geo.lum_disk_back = qdisk.heat[kkk] += pp.w; // 60a - ksl - Added to be able to calculate illum of disk + qdisk.ave_freq[kkk] += pp.w * pp.freq; + + if (geo.absorb_reflect == BACK_RAD_SCATTER) + { + randvcos (pp.lmn, normal); + stuff_phot (&pp, p); + tau_scat = -log (1. - (rand () + 0.5) / MAXRAND); + istat = pp.istat = P_INWIND; // if we got here, the photon stays in the wind- make sure istat doesn't say scattered still! + tau = 0; + if (iextract) + { + stuff_phot (&pp, &pextract); + extract (w, &pextract, PTYPE_DISK); // Treat as wind photon for purpose of extraction + } + } + else + { /*This is the end of the line for this photon */ + stuff_phot (&pp, p); + break; + } + } + + if (istat == P_SCAT) + { /* Cause the photon to scatter and reinitilize */ + + /* 71 - 1112 - ksl - placed this line here to try to avoid an error I was seeing in scatter. I believe the first if + statement has a loophole that needs to be plugged, when it comes back with avalue of n = -1 */ + + pp.grid = n = where_in_grid (wmain[pp.grid].ndom, pp.x); + + if (n < 0) + { + Error + ("trans_phot: Trying to scatter a photon which is not in the wind\n"); + Error ("trans_phot: grid %3d x %8.2e %8.2e %8.2e\n", pp.grid, + pp.x[0], pp.x[1], pp.x[2]); + Error ("trans_phot: This photon is effectively lost!\n"); + istat = pp.istat = p->istat = P_ERROR; + stuff_phot (&pp, p); + break; + } + + /* 57+ -- ksl -- Add check to see if there is a cell in the plasma structure for this. This is a problem that needs + fixing */ + /* 1506 JM -- this appeared to happen due to a rather convoluted problem involving DFUDGE + and not updating the istat variable properly. See Issue #154 for discussion */ + + if (wmain[n].nplasma == NPLASMA) + { + Error + ("trans_phot: Trying to scatter a photon which is not in a cell in the plasma structure\n"); + Error ("trans_phot: grid %3d x %8.2e %8.2e %8.2e\n", pp.grid, + pp.x[0], pp.x[1], pp.x[2]); + Error ("trans_phot: This photon is effectively lost!\n"); + istat = pp.istat = p->istat = P_ERROR; + stuff_phot (&pp, p); + break; + } + + /* 57h -- ksl -- Add check to verify the cell has some volume */ + + if (wmain[n].vol <= 0) + { + Error + ("trans_phot: Trying to scatter a photon in a cell with no wind volume\n"); + Error ("trans_phot: grid %3d x %8.2e %8.2e %8.2e\n", pp.grid, + pp.x[0], pp.x[1], pp.x[2]); + Log ("istat %d\n", pp.istat); + Error ("trans_phot: This photon is effectively lost!\n"); + istat = pp.istat = p->istat = P_ERROR; + stuff_phot (&pp, p); + break; + + } + + /* 0215 SWM - Added cell-based reverberation mapping */ + if ((geo.reverb == REV_WIND || geo.reverb == REV_MATOM) + && geo.ioniz_or_extract && geo.wcycle == geo.wcycles - 1) + { + wind_paths_add_phot (&wmain[n], &pp); + } + + + /* SS July 04 - next lines modified so that the "thermal trapping" model of anisotropic scattering is included in the + macro atom method. What happens now is all in scatter - within that routine the "thermal trapping" model is used to + decide what the direction of emission is before returning here. 54b-ksl -- To see what the code did previously see + py46. I've confirmed that the current version of scattering really does what the old code did for two-level lines */ + + + nnscat = 0; + nnscat++; + ptr_nres = &nres; + + /* 74a_ksl - Check added to search for error in weights */ + if (sane_check (pp.w)) + { + Error + ("trans_phot:sane_check photon %d has weight %e before scatter\n", + p->np, pp.w); + } + if ((nerr = scatter (&pp, ptr_nres, &nnscat)) != 0) + { + Error ("trans_phot: Bad return from scatter %d at point 2", + nerr); + } + pp.nscat++; + /* 74a_ksl - Check added to search for error in weights */ + + if (sane_check (pp.w)) + { + Error + ("trans_phot:sane_check photon %d has weight %e after scatter\n", + p->np, pp.w); + } + + /* SS June 04: During the spectrum calculation cycles, photons are thrown away when they interact with macro atoms or + become k-packets. This is done by setting their weight to zero (effectively they have been absorbed into either + excitation energy or thermal energy). Since they now have no weight there is no need to follow them further. */ + /* 54b-ksl ??? Stuart do you really mean the comment above; it's not obvious to me since if true why does one need to + calculate the progression of photons through the wind at all??? Also how is this enforced; where is pp.w set to a + low value. */ + /* JM 1504 -- This is correct. It's one of the odd things about combining the macro-atom approach with our way of doing + 'spectral cycles'. If photons activate macro-atoms they are destroyed, but we counter this by generating photons + from deactivating macro-atoms with the already calculated emissivities. */ + + if (geo.matom_radiation == 1 && geo.rt_mode == RT_MODE_MACRO + && pp.w < weight_min) + /* Flag for the spectrum calculations in a macro atom calculation SS */ + { + istat = pp.istat = P_ABSORB; + pp.tau = VERY_BIG; + stuff_phot (&pp, p); + break; + } + + // Calculate the line heating and if the photon was absorbed break finish up + // ??? Need to modify line_heat for multiple scattering but not yet + // Condition that nres < nlines added (SS) + + if (nres > -1 && nres < nlines) + { + pp.nrscat++; + + /* This next statement writes out the position of every resonant scattering event to a file */ + if (modes.track_resonant_scatters) + fprintf (pltptr, + "Photon %i has resonant scatter at %.2e %.2e %.2e in wind cell %i (grid cell=%i). Freq=%e Weight=%e\n", + p->np, pp.x[0], pp.x[1], pp.x[2], wmain[n].nplasma, + pp.grid, pp.freq, pp.w); + + /* 68a - 090124 - ksl - Increment the number of scatters by this ion in this cell */ + /* 68c - 090408 - ksl - Changed this to the weight of the photon at the time of the scatter */ + + plasmamain[wmain[n].nplasma].scatters[line[nres].nion] += pp.w; + + if (geo.rt_mode == RT_MODE_2LEVEL) // only do next line for non-macro atom case + { + line_heat (&plasmamain[wmain[n].nplasma], &pp, nres); + } + + if (pp.w < weight_min) + { + istat = pp.istat = P_ABSORB; /* This photon was absorbed by continuum opacity within the wind */ + pp.tau = VERY_BIG; + stuff_phot (&pp, p); + break; + } + } + + + /* The next if statement causes photons to be extracted during the creation of the detailed spectrum portion of the + program */ + + /* N.B. To use the anisotropic scattering option, extract needs to follow scatter. This is because the reweighting + which occurs in extract needs the pdf for scattering to have been initialized. 02may ksl. This seems to be OK at + present. */ + + if (iextract) + { + stuff_phot (&pp, &pextract); + + + /* JM 1407 -- This next loop is required because in anisotropic scattering mode 2 we have normalised our rejection + method. This means that we have to adjust nnscat by this factor, since nnscat will be lower by a factor of + 1/p_norm */ + if (geo.scatter_mode == SCATTER_MODE_THERMAL + && pextract.nres <= NLINES && pextract.nres > 0) + { + /* we normalised our rejection method by the escape probability along the vector of maximum velocity gradient. + First find the sobolev optical depth along that vector */ + tau_norm = + sobolev (&wmain[pextract.grid], pextract.x, -1.0, + lin_ptr[pextract.nres], + wmain[pextract.grid].dvds_max); + + /* then turn into a probability */ + p_norm = p_escape_from_tau (tau_norm); + + } + else + { + p_norm = 1.0; + + /* throw an error if nnscat does not equal 1 */ + /* JM 1504-- originally I'd used the wrong value of nnscat here. This would throw large amounts of errors which + weren't actually errors */ + /* nnscat is the quantity associated with this photon being extracted */ + if (nnscat != 1) + Error + ("nnscat is %i for photon %i in scatter mode %i! nres %i NLINES %i\n", + nnscat, p->np, geo.scatter_mode, pextract.nres, + NLINES); + } + + /* We then increase weight to account for number of scatters. This is done because in extract we multiply by the + escape probability along a given direction, but we also need to divide the weight by the mean escape + probability, which is equal to 1/nnscat */ + pextract.w *= nnscat / p_norm; + + if (sane_check (pextract.w)) + { + Error + ("trans_phot: sane_check photon %d has weight %e before extract\n", + p->np, pextract.w); + } + extract (w, &pextract, PTYPE_WIND); // Treat as wind photon for purpose of extraction + } + + + + + /* OK we are ready to continue the processing of a photon which has scattered. + * The next steps reinitialize parameters + so that the photon can continue throug the wind */ + + tau_scat = -log (1. - (rand () + 0.5) / MAXRAND); + istat = pp.istat = P_INWIND; // if we got here, the photon stays in the wind- make sure istat doesn't say scattered still! + tau = 0; + + stuff_v (pp.x, x_dfudge_check); // this is a vector we use to see if dfudge moved the photon outside the wind cone + reposition (&pp); + + /* JM 1506 -- call walls again to account for instance where DFUDGE + can take photon outside of the wind and into the disk or star + after scattering. Note that walls updates the istat in pp as well. + This may not be necessary but I think to account for every eventuality + it should be done */ + istat = walls (&pp, p, normal); + + /* This *does not* update istat if the photon scatters outside of the wind- + I guess P_INWIND is really in wind or empty space but not escaped. + translate_in_space will take care of this next time round. All a bit + convoluted but should work. */ + + /* JM 1506 -- we don't throw errors here now, but we do keep a track + of how many 4 photons were lost due to DFUDGE pushing them + outside of the wind after scatter */ + + // XXX PLACEHOLDER Check that this is the correct logic here + if (where_in_wind (pp.x, &ndom) != W_ALL_INWIND + && where_in_wind (x_dfudge_check, &ndom) == W_ALL_INWIND) + { + n_lost_to_dfudge++; // increment the counter (checked at end of trans_phot) + } + + stuff_phot (&pp, p); + icell = 0; + } + + + + /* This completes the portion of the code that handles the scattering of a photon. + * What follows is a simple check to see if + this particular photon has gotten stuck in the wind 54b-ksl */ + + if (pp.nscat == MAXSCAT) + { + istat = pp.istat = P_TOO_MANY_SCATTERS; /* Scattered too many times */ + stuff_phot (&pp, p); + break; + } + + if (pp.istat == P_ADIABATIC) + { + istat = pp.istat = p->istat = P_ADIABATIC; + stuff_phot (&pp, p); + break; + } + + /* This appears partly to be an insurance policy. It is not obvious that for example nscat + * and nrscat need to be updated */ + + p->istat = istat; + p->nscat = pp.nscat; + p->nrscat = pp.nrscat; + p->w = pp.w; // Assure that final weight of photon is returned. - if (istat == P_HIT_STAR) - { /* It hit the star */ - geo.lum_star_back+=pp.w; - if (geo.absorb_reflect==BACK_RAD_SCATTER){ - randvcos(pp.lmn,normal); - stuff_phot (&pp, p); - tau_scat = -log (1. - (rand () + 0.5) / MAXRAND); - istat = pp.istat = P_INWIND; // if we got here, the photon stays in the wind- make sure istat doesn't say scattered still! - tau = 0; - if (iextract) { - stuff_phot (&pp, &pextract); - extract (w, &pextract, PTYPE_STAR); // Treat as wind photon for purpose of extraction - } - } - else { /*This is the end of the line for this photon */ - stuff_phot (&pp, p); - break; - } } + /* This is the end of the loop over a photon */ - if (istat == P_HIT_DISK) - { - /* It hit the disk */ - - /* Store the energy of the photon bundle into a disk structure so that one can determine later how much and where the - disk was heated by photons. - Note that the disk is defined from 0 to NRINGS-2. NRINGS-1 contains the position of the outer radius of the disk. */ - - rrr = sqrt (dot (pp.x, pp.x)); - kkk = 0; - while (rrr > qdisk.r[kkk] && kkk < NRINGS - 1) - kkk++; - kkk--; // So that the heating refers to the heating between kkk and kkk+1 - qdisk.nhit[kkk]++; - geo.lum_disk_back=qdisk.heat[kkk] += pp.w; // 60a - ksl - Added to be able to calculate illum of disk - qdisk.ave_freq[kkk] += pp.w * pp.freq; - - if (geo.absorb_reflect==BACK_RAD_SCATTER){ - randvcos(pp.lmn,normal); - stuff_phot (&pp, p); - tau_scat = -log (1. - (rand () + 0.5) / MAXRAND); - istat = pp.istat = P_INWIND; // if we got here, the photon stays in the wind- make sure istat doesn't say scattered still! - tau = 0; - if (iextract) { - stuff_phot (&pp, &pextract); - extract (w, &pextract, PTYPE_DISK); // Treat as wind photon for purpose of extraction - } - } - else { /*This is the end of the line for this photon */ - stuff_phot (&pp, p); - break; - } - } + /* The next section is for diagnostic purpposes. There are two possibilities. If you wish to know where + * the photon was last while in the wind, you want to track p; if you wish to know where it hits the + * outer boundary of the calculation you would want pp */ - if (istat == P_SCAT) - { /* Cause the photon to scatter and reinitilize */ - - /* 71 - 1112 - ksl - placed this line here to try to avoid an error I was seeing in scatter. I believe the first if - statement has a loophole that needs to be plugged, when it comes back with avalue of n = -1 */ - - pp.grid = n = where_in_grid (wmain[pp.grid].ndom, pp.x); - - if (n < 0) - { - Error ("trans_phot: Trying to scatter a photon which is not in the wind\n"); - Error ("trans_phot: grid %3d x %8.2e %8.2e %8.2e\n", pp.grid, pp.x[0], pp.x[1], pp.x[2]); - Error ("trans_phot: This photon is effectively lost!\n"); - istat = pp.istat = p->istat = P_ERROR; - stuff_phot (&pp, p); - break; - } - - /* 57+ -- ksl -- Add check to see if there is a cell in the plasma structure for this. This is a problem that needs - fixing */ - /* 1506 JM -- this appeared to happen due to a rather convoluted problem involving DFUDGE - and not updating the istat variable properly. See Issue #154 for discussion */ - - if (wmain[n].nplasma == NPLASMA) - { - Error ("trans_phot: Trying to scatter a photon which is not in a cell in the plasma structure\n"); - Error ("trans_phot: grid %3d x %8.2e %8.2e %8.2e\n", pp.grid, pp.x[0], pp.x[1], pp.x[2]); - Error ("trans_phot: This photon is effectively lost!\n"); - istat = pp.istat = p->istat = P_ERROR; - stuff_phot (&pp, p); - break; - } - - /* 57h -- ksl -- Add check to verify the cell has some volume */ - - if (wmain[n].vol <= 0) - { - Error ("trans_phot: Trying to scatter a photon in a cell with no wind volume\n"); - Error ("trans_phot: grid %3d x %8.2e %8.2e %8.2e\n", pp.grid, pp.x[0], pp.x[1], pp.x[2]); - Log ("istat %d\n", pp.istat); - Error ("trans_phot: This photon is effectively lost!\n"); - istat = pp.istat = p->istat = P_ERROR; - stuff_phot (&pp, p); - break; - - } - - /* 0215 SWM - Added cell-based reverberation mapping */ - if ((geo.reverb == REV_WIND || geo.reverb == REV_MATOM) && geo.ioniz_or_extract && geo.wcycle == geo.wcycles - 1) - { - wind_paths_add_phot (&wmain[n], &pp); - } - - - /* SS July 04 - next lines modified so that the "thermal trapping" model of anisotropic scattering is included in the - macro atom method. What happens now is all in scatter - within that routine the "thermal trapping" model is used to - decide what the direction of emission is before returning here. 54b-ksl -- To see what the code did previously see - py46. I've confirmed that the current version of scattering really does what the old code did for two-level lines */ - - - nnscat = 0; - nnscat++; - ptr_nres = &nres; - - /* 74a_ksl - Check added to search for error in weights */ - if (sane_check (pp.w)) - { - Error ("trans_phot:sane_check photon %d has weight %e before scatter\n", p->np, pp.w); - } - if ((nerr = scatter (&pp, ptr_nres, &nnscat)) != 0) - { - Error ("trans_phot: Bad return from scatter %d at point 2", nerr); - } - pp.nscat++; - /* 74a_ksl - Check added to search for error in weights */ - - if (sane_check (pp.w)) - { - Error ("trans_phot:sane_check photon %d has weight %e after scatter\n", p->np, pp.w); - } - - /* SS June 04: During the spectrum calculation cycles, photons are thrown away when they interact with macro atoms or - become k-packets. This is done by setting their weight to zero (effectively they have been absorbed into either - excitation energy or thermal energy). Since they now have no weight there is no need to follow them further. */ - /* 54b-ksl ??? Stuart do you really mean the comment above; it's not obvious to me since if true why does one need to - calculate the progression of photons through the wind at all??? Also how is this enforced; where is pp.w set to a - low value. */ - /* JM 1504 -- This is correct. It's one of the odd things about combining the macro-atom approach with our way of doing - 'spectral cycles'. If photons activate macro-atoms they are destroyed, but we counter this by generating photons - from deactivating macro-atoms with the already calculated emissivities. */ - - if (geo.matom_radiation == 1 && geo.rt_mode == RT_MODE_MACRO && pp.w < weight_min) - /* Flag for the spectrum calculations in a macro atom calculation SS */ - { - istat = pp.istat = P_ABSORB; - pp.tau = VERY_BIG; - stuff_phot (&pp, p); - break; - } - - // Calculate the line heating and if the photon was absorbed break finish up - // ??? Need to modify line_heat for multiple scattering but not yet - // Condition that nres < nlines added (SS) - - if (nres > -1 && nres < nlines) - { - pp.nrscat++; - - /* This next statement writes out the position of every resonant scattering event to a file */ - if (modes.track_resonant_scatters) - fprintf (pltptr, - "Photon %i has resonant scatter at %.2e %.2e %.2e in wind cell %i (grid cell=%i). Freq=%e Weight=%e\n", - p->np, pp.x[0], pp.x[1], pp.x[2], wmain[n].nplasma, pp.grid, pp.freq, pp.w); - - /* 68a - 090124 - ksl - Increment the number of scatters by this ion in this cell */ - /* 68c - 090408 - ksl - Changed this to the weight of the photon at the time of the scatter */ - - plasmamain[wmain[n].nplasma].scatters[line[nres].nion] += pp.w; - - if (geo.rt_mode == RT_MODE_2LEVEL) // only do next line for non-macro atom case - { - line_heat (&plasmamain[wmain[n].nplasma], &pp, nres); - } - - if (pp.w < weight_min) - { - istat = pp.istat = P_ABSORB; /* This photon was absorbed by continuum opacity within the wind */ - pp.tau = VERY_BIG; - stuff_phot (&pp, p); - break; - } - } - - - /* The next if statement causes photons to be extracted during the creation of the detailed spectrum portion of the - program */ - - /* N.B. To use the anisotropic scattering option, extract needs to follow scatter. This is because the reweighting - which occurs in extract needs the pdf for scattering to have been initialized. 02may ksl. This seems to be OK at - present. */ - - if (iextract) - { - stuff_phot (&pp, &pextract); - - - /* JM 1407 -- This next loop is required because in anisotropic scattering mode 2 we have normalised our rejection - method. This means that we have to adjust nnscat by this factor, since nnscat will be lower by a factor of - 1/p_norm */ - if (geo.scatter_mode == SCATTER_MODE_THERMAL && pextract.nres <= NLINES && pextract.nres > 0) - { - /* we normalised our rejection method by the escape probability along the vector of maximum velocity gradient. - First find the sobolev optical depth along that vector */ - tau_norm = sobolev (&wmain[pextract.grid], pextract.x, -1.0, lin_ptr[pextract.nres], wmain[pextract.grid].dvds_max); - - /* then turn into a probability */ - p_norm = p_escape_from_tau (tau_norm); - - } - else - { - p_norm = 1.0; - - /* throw an error if nnscat does not equal 1 */ - /* JM 1504-- originally I'd used the wrong value of nnscat here. This would throw large amounts of errors which - weren't actually errors */ - /* nnscat is the quantity associated with this photon being extracted */ - if (nnscat != 1) - Error - ("nnscat is %i for photon %i in scatter mode %i! nres %i NLINES %i\n", - nnscat, p->np, geo.scatter_mode, pextract.nres, NLINES); - } - - /* We then increase weight to account for number of scatters. This is done because in extract we multiply by the - escape probability along a given direction, but we also need to divide the weight by the mean escape - probability, which is equal to 1/nnscat */ - pextract.w *= nnscat / p_norm; - - if (sane_check (pextract.w)) - { - Error ("trans_phot: sane_check photon %d has weight %e before extract\n", p->np, pextract.w); - } - extract (w, &pextract, PTYPE_WIND); // Treat as wind photon for purpose of extraction - } - - - - - /* OK we are ready to continue the processing of a photon which has scattered. The next steps reinitialize parameters - so that the photon can continue throug the wind */ - - tau_scat = -log (1. - (rand () + 0.5) / MAXRAND); - istat = pp.istat = P_INWIND; // if we got here, the photon stays in the wind- make sure istat doesn't say scattered still! - tau = 0; - - stuff_v (pp.x, x_dfudge_check); // this is a vector we use to see if dfudge moved the photon outside the wind cone - reposition (&pp); - - /* JM 1506 -- call walls again to account for instance where DFUDGE - can take photon outside of the wind and into the disk or star - after scattering. Note that walls updates the istat in pp as well. - This may not be necessary but I think to account for every eventuality - it should be done */ - istat = walls (&pp, p,normal); - - /* This *does not* update istat if the photon scatters outside of the wind- - I guess P_INWIND is really in wind or empty space but not escaped. - translate_in_space will take care of this next time round. All a bit - convoluted but should work. */ - - /* JM 1506 -- we don't throw errors here now, but we do keep a track - of how many 4 photons were lost due to DFUDGE pushing them - outside of the wind after scatter */ - - // XXX PLACEHOLDER Check that this is the correct logic here - if (where_in_wind (pp.x, &ndom) != W_ALL_INWIND && where_in_wind (x_dfudge_check, &ndom) == W_ALL_INWIND) - { - n_lost_to_dfudge++; // increment the counter (checked at end of trans_phot) - } - - stuff_phot (&pp, p); - icell = 0; - } - - - - /* This completes the portion of the code that handles the scattering of a photon What follows is a simple check to see if - this particular photon has gotten stuck in the wind 54b-ksl */ - - if (pp.nscat == MAXSCAT) + if (modes.save_photons) { - istat = pp.istat = P_TOO_MANY_SCATTERS; /* Scattered too many times */ - stuff_phot (&pp, p); - break; + // save_photons (p, "Final"); // Where the last position of the photon in the wind + save_photons (&pp, "Final"); //The postion of the photon where it exits the calculation } - - if (pp.istat == P_ADIABATIC) - { - istat = pp.istat = p->istat = P_ADIABATIC; - stuff_phot (&pp, p); - break; - } - - /* This appears partly to be an insurance policy. It is not obvious that for example nscat and nrscat need to be updated */ - p->istat = istat; - p->nscat = pp.nscat; - p->nrscat = pp.nrscat; - p->w = pp.w; // Assure that final weight of photon is returned. - - } - /* This is the end of the loop over individual photons */ return (0); }