Merge pull request #167 from Higginbottom/dev

Compton Scattering

source/compton.c

@@ -27,6 +27,8 @@
 #include <stdlib.h>
 #include "atomic.h"
 #include "python.h"
+#include <gsl/gsl_rng.h>
+
 
 /**************************************************************************
                     Space Telescope Science Institute
@@ -73,6 +75,8 @@ kappa_comp (xplasma, freq)
   return (x);
 }
 
+
+
 /**************************************************************************
                     Space Telescope Science Institute
 
@@ -230,3 +234,193 @@ klein_nishina (nu)
 
   return (kn);
 }
+
+double z_rand,sigma_tot,x1; //External variables to allow zfunc to search for the correct fractional energy change
+
+/**************************************************************************
+                    Southampton University
+
+
+  Synopsis:  compton_dir computes a random direction (and hence frequency change)
+			for a photon undergoing compton scattering. At low frequencies, this
+			is just thompson scattering.
+
+  Description:	
+
+  Arguments:  
+			p - the photon currently being scattered, this gives us the current direction and the frequency
+			xplasma- pointer to the current plasma cell
+
+  Returns:   kappa - nothing, but updates the photon frequency and direction
+
+  Notes:   
+
+  History:
+2015	NSH coded as part of teh summer 2015 code sprint
+
+ ************************************************************************/
+
+
+
+int
+	compton_dir (p,xplasma)
+		PhotPtr p;
+        PlasmaPtr xplasma;		// Pointer to current plasma cell
+
+{
+	double f_min,f_max,f; //The theoretical maxmimum energy change
+	double n,l,m,phi,len;  //The direction cosines of the new photon direction in the frame of reference with q along the photon path
+	struct basis nbasis;  //The basis function which transforms between the photon frame and the pbserver frame	
+	double lmn[3]; /* the individual direction cosines in the rotated frame */
+	double x[3]; /*photon direction in the frame of reference of the original photon*/
+	double dummy[3],c[3];
+
+	x1=H*p->freq/MELEC/C/C; //compute the ratio of photon energy to electron energy. In the electron rest frame this is just the electron rest mass energ 
+
+	n=l=m=0.0;  //initialise some variables to avoid warnings
+
+
+	if (x1<0.0001) //If the photon energy is much less than electron mass, we just have thompson scattering.
+	{
+		randvec(lmn,1.0); //Generate a normal isotropic scatter
+		f=1.0;   //There is no energy loss
+		stuff_v (lmn, p->lmn);
+	}
+	else
+	{
+		z_rand=rand()/MAXRAND; //Generate a random number between 0 and 1 - this is the random location in the klein nishina scattering distribution - it gives the energy loss and also direction.
+		f_min=1.;         //The minimum energy loss - i.e. no energy loss
+		f_max=1.+(2.*x1);  //The maximum energy loss
+
+		sigma_tot=sigma_compton_partial(f_max,x1);  //Communicated externally to the integrand function in the zbrent call below, this is the maximum cross section, used to scale the K_N function to lie between 0 and 1.
+
+		f=zbrent(compton_func,f_min,f_max,1e-8);  //Find the zero point of the function compton_func - this finds the point in the KN function that represents our random energy loss.
+		n=(1.-((f-1.)/x1));   //This is the angle cosine of the new direction in the frame of reference of the photon
+//		printf ("f=%e n=%e fmin=%e fmax=%e\n",f,n,f_min,f_max);
+
+		if (isfinite(len=sqrt(1.-(n*n)))==0)   //Compute the length of the other sides - the isfinite is to take care of the very rare occasion where n=1!
+		len=0.0;
+		phi = 0.0;           //no need to randomise phi, the random rotation of the vector generating the basis function takes care of this
+
+		l = len*cos (phi);   //compute the angle cosines of the other two dimensions.
+		m = len*sin (phi);
+
+		randvec(dummy,1.0);  //Get a random vector
+
+		cross(dummy, p->lmn, c);  //c will be perpendicular to p->lmn
+
+//		printf ("c= %e %e %e n=%e acos(n)=%f\n",c[0],c[1],c[2],n,acos(n)*RADIAN);
+
+		create_basis (p->lmn, c, &nbasis); //create a basis with the first axis in the direction of the original photon direction, c will be perpendicular to the photon direction. Orientaion of the y/z axes are will give the randomization of the phi axis
+
+
+		x[0] = n;  //This is the cosine direction of the new photon direction, in the frame of reference where the first axis is in the original photon direction
+		x[1] = l;
+		x[2] = m;
+
+
+      project_from (&nbasis, x, lmn);	/* Project the vector from the FOR of the original photon into the observer frame */
+		renorm(lmn,1.0);
+//		Log("f=%e freq=%e  n=%e l_old=%e m_old=%e n_old=%e l_new=%e m_new=%e n_new=%e len=%e\n",f,p->freq,n,p->lmn[0],p->lmn[1],p->lmn[2],lmn[0],lmn[1],lmn[2],length(lmn));
+//	  Log_flush();
+//	  renorm (lmn,1.0);
+
+		stuff_v (lmn,p->lmn);
+//		randvec(a,1.0); //Generate a normal isotropic scatter
+//				f=1.0;   //There is no energy loss
+//		stuff_v (a, p->lmn);
+//	printf ("Original %e %e %e theta %e new %e %e %e\n",pold.lmn[0],pold.lmn[1],pold.lmn[2],n,p->lmn[0],p->lmn[1],p->lmn[2]);		
+
+	}
+
+	p->freq=p->freq/f;  //reduce the photon frequency
+	p->w=p->w/f;  //reduce the photon weight by the same ammount to conserve photon numbers
+return(0);
+}
+
+/**************************************************************************
+                    Southampton University
+
+
+  Synopsis:  compton_func is a simple function that is equal to zero when the 
+			external variable z_rand is equal to the normalised KN cross section.
+			It is used in a call to zbrent in the function compton_dir
+
+  Description:	
+
+  Arguments:  
+			f- the fractional energy change which is supplied to sigma_compton_partial
+
+  Returns:   the difference between sigma(f) and the randomised cross section we are searching for
+
+  Notes:   Since this function is called by zbrent, many of the variables have to be communicated
+			externally. These are 
+				x1, the ratio of photon energy to electron rest mass,
+				sigma_tot, the total angle integrated KN cross section
+				z_rand - a randomised number between 0 and 1 representing the normalised cross section we want to find
+		
+
+  History:
+2015	NSH coded as part of the summer 2015 code sprint
+
+ ************************************************************************/
+
+
+double
+	compton_func(f)
+		double f;
+{
+	double ans;
+	ans=(sigma_compton_partial(f,x1)/sigma_tot)-z_rand;
+		return(ans);
+}
+
+/**************************************************************************
+                    Southampton University
+
+
+  Synopsis:  sigma_compton_partial is the KN cross section as a function of the 
+			fractional energy change
+
+  Description:	
+
+  Arguments:  
+			f - the fractional energy change
+			x	the enerrgy of the photon divided by the rest mass energy of an electron
+
+  Returns:   the cross section for the scattering angle that gives this fractional energy change.
+
+  Notes:   
+			Stuart Sim supplied this formula, but coundn't recall where he had found it, although
+			he had rederived it and thinks it is correct!
+
+  History:
+2015	NSH coded as part of the summer 2015 code sprint
+
+ ************************************************************************/
+
+
+double 
+   sigma_compton_partial(f,x)
+      double f; //This is the fractional energy change, nu/nu'
+   double x; //h nu/mec**2 - the energy of the photon divided by the rest energy of an eectron
+   {
+      double term1,term2,term3,tot;
+
+   term1 = ( (x*x) - (2*x) - 2 ) * log(f) / x / x;
+   term2 = ( ((f*f) -1) / (f * f)) / 2;
+   term3 = ( (f - 1) / x) * ( (1/x) + (2/f) + (1/(x*f)));
+
+   tot = 3 * THOMPSON * (term1 + term2 + term3) / (8 * x);
+
+   return(tot);
+
+}
+
+
+
+
+
+
+
+

source/foo.h

@@ -422,6 +422,9 @@ double kappa_comp(PlasmaPtr xplasma, double freq);
 double kappa_ind_comp(PlasmaPtr xplasma, double freq);
 double total_comp(WindPtr one, double t_e);
 double klein_nishina(double nu);
+int compton_dir(PhotPtr p, PlasmaPtr xplasma);
+double compton_func(double f);
+double sigma_compton_partial(double f, double x);
 /* torus.c */
 double torus_rho(double x[]);
 double ds_to_cylinder(double rho, struct photon *p);

source/knigge.c

@@ -87,7 +87,7 @@ Terminolgy is awful here. -- ksl
 As now represented geo.kn_dratio is the distance to the focus point in stellar radii!
 
 */
-  rddoub ("kn.d", &geo.kn_dratio);
+  rddoub ("kn.d(in_wd_radii)", &geo.kn_dratio);
   Log_silent ("dmin = %f so the ratio d/dmin here is %f  (%.2e %.2e) \n",
 	      dmin, geo.kn_dratio / dmin, geo.diskrad, geo.rstar);
 
@@ -112,8 +112,8 @@ to be modified -- ksl 04jun */
   /* JM 1502 -- added capability for user to adjust launch radii of KWD wind
      required to reproduce Stuart's X-ray models (Sim+ 2008,2010) 
      units are in stellar radii / WD radii / grav radii as in SV model */
-  rddoub ("kn.rmin", &geo.wind_rmin);
-  rddoub ("kn.rmax", &geo.wind_rmax ); 
+  rddoub ("kn.rmin(in_wd_radii)", &geo.wind_rmin);
+  rddoub ("kn.rmax(in_wd_radii)", &geo.wind_rmax ); 
 
   geo.wind_rmin *= geo.rstar;
   geo.wind_rmax *= geo.rstar;

source/para_update.c

@@ -48,7 +48,7 @@ communicate_estimators_para ()
 
   /* The size of the helper array for doubles. We transmit 10 numbers 
      for each cell, plus three arrays, each of length NXBANDS */
-  plasma_double_helpers = (14 + 3 * NXBANDS) * NPLASMA;
+  plasma_double_helpers = (15 + 3 * NXBANDS) * NPLASMA;
 
   /* The size of the helper array for integers. We transmit 6 numbers 
      for each cell, plus one array of length NXBANDS */
@@ -89,11 +89,13 @@ communicate_estimators_para ()
       redhelper[mpi_i + 11 * NPLASMA] = plasmamain[mpi_i].j_scatt / np_mpi_global;
       redhelper[mpi_i + 12 * NPLASMA] = plasmamain[mpi_i].ip_direct / np_mpi_global;
       redhelper[mpi_i + 13 * NPLASMA] = plasmamain[mpi_i].ip_scatt / np_mpi_global;
+      redhelper[mpi_i + 14 * NPLASMA] = plasmamain[mpi_i].heat_auger / np_mpi_global;
+
       for (mpi_j = 0; mpi_j < NXBANDS; mpi_j++)
 	{
-	  redhelper[mpi_i + (14 + mpi_j) * NPLASMA] = plasmamain[mpi_i].xj[mpi_j] / np_mpi_global;
-	  redhelper[mpi_i + (14 + NXBANDS + mpi_j) * NPLASMA] = plasmamain[mpi_i].xave_freq[mpi_j] / np_mpi_global;
-	  redhelper[mpi_i + (14 + 2 * NXBANDS + mpi_j) * NPLASMA] = plasmamain[mpi_i].xsd_freq[mpi_j] / np_mpi_global;
+	  redhelper[mpi_i + (15 + mpi_j) * NPLASMA] = plasmamain[mpi_i].xj[mpi_j] / np_mpi_global;
+	  redhelper[mpi_i + (15 + NXBANDS + mpi_j) * NPLASMA] = plasmamain[mpi_i].xave_freq[mpi_j] / np_mpi_global;
+	  redhelper[mpi_i + (15 + 2 * NXBANDS + mpi_j) * NPLASMA] = plasmamain[mpi_i].xsd_freq[mpi_j] / np_mpi_global;
 
 	  /* 131213 NSH populate the band limited min and max frequency arrays */
 	  maxbandfreqhelper[mpi_i * NXBANDS + mpi_j] = plasmamain[mpi_i].fmax[mpi_j];
@@ -136,11 +138,13 @@ communicate_estimators_para ()
       plasmamain[mpi_i].j_scatt = redhelper2[mpi_i + 11 * NPLASMA];
       plasmamain[mpi_i].ip_direct = redhelper2[mpi_i + 12 * NPLASMA];
       plasmamain[mpi_i].ip_scatt = redhelper2[mpi_i + 13 * NPLASMA];
+      plasmamain[mpi_i].heat_auger = redhelper2[mpi_i + 14 * NPLASMA];
+
       for (mpi_j = 0; mpi_j < NXBANDS; mpi_j++)
 	{
-	  plasmamain[mpi_i].xj[mpi_j] = redhelper2[mpi_i + (14 + mpi_j) * NPLASMA];
-	  plasmamain[mpi_i].xave_freq[mpi_j] = redhelper2[mpi_i + (14 + NXBANDS + mpi_j) * NPLASMA];
-	  plasmamain[mpi_i].xsd_freq[mpi_j] = redhelper2[mpi_i + (14 + NXBANDS * 2 + mpi_j) * NPLASMA];
+	  plasmamain[mpi_i].xj[mpi_j] = redhelper2[mpi_i + (15 + mpi_j) * NPLASMA];
+	  plasmamain[mpi_i].xave_freq[mpi_j] = redhelper2[mpi_i + (15 + NXBANDS + mpi_j) * NPLASMA];
+	  plasmamain[mpi_i].xsd_freq[mpi_j] = redhelper2[mpi_i + (15 + NXBANDS * 2 + mpi_j) * NPLASMA];
 
 	  /* 131213 NSH And unpack the min and max banded frequencies to the plasma array */
 	  plasmamain[mpi_i].fmax[mpi_j] = maxbandfreqhelper2[mpi_i * NXBANDS + mpi_j];

source/photon2d.c

@@ -58,7 +58,10 @@ translate (w, pp, tau_scat, tau, nres)
     }
   else if ((pp->grid = where_in_grid (pp->x)) >= 0)
     {
+//		 printf ("photon %i start=%e %e %e %e",pp->np,pp->x[0], pp->x[1], pp->x[2],sqrt(pp->x[0]*pp->x[0]+pp->x[1]*pp->x[1]+pp->x[2]*pp->x[2]));
       istat = translate_in_wind (w, pp, tau_scat, tau, nres);
+//	 printf ("end=%e %e %e %e\n",pp->x[0], pp->x[1], pp->x[2],sqrt(pp->x[0]*pp->x[0]+pp->x[1]*pp->x[1]+pp->x[2]*pp->x[2]));
+
     }
   else
     {

source/radiation.c

@@ -96,6 +96,8 @@
 	1405    JM corrected error (#73) so that photon frequency is shifted to the rest frame of the 
 	        cell in question. Also added check which checks if a photoionization edge is crossed
 	        along ds.
+	1508	NSH slight modification to mean that compton scattering no longer reduces the weight of
+			the photon in this part of the code. It is now done when the photon scatters.
 **************************************************************/
 
 #include <stdio.h>
@@ -367,24 +369,53 @@ if (freq > phot_freq_min)
       Error ("Radiation:sane_check CHECKING ff=%e, comp=%e, ind_comp=%e\n",
 	     frac_ff, frac_comp, frac_ind_comp);
     }
-/* Calculate the reduction in the w of the photon bundle along with the average
-   weight in the cell */
+/* Calculate the heating effect*/
 
   if (tau > 0.0001)
     {				/* Need differentiate between thick and thin cases */
       x = exp (-tau);
-      p->w = w_out = w_in * x;
-      energy_abs = w_in - w_out;
-      w_ave = (w_in - w_out) / tau;
+      energy_abs = w_in *(1.-x);
     }
   else
     {
       tau2 = tau * tau;
-      p->w = w_out = w_in * (1. - tau + 0.5 * tau2);	/*Calculate to second order */
       energy_abs = w_in * (tau - 0.5 * tau2);
-      w_ave = w_in * (1. - 0.5 * tau + 0.1666667 * tau2);
     }
 
+	/* Calculate the reduction in weight - compton scattering is not included, it is now included at scattering
+	however induced compton heating is not implemented at scattering, so it should remain here for the time being
+	to maimtain consistency.*/
+
+    tau = (kappa_tot-frac_comp) * ds;
+
+    if (sane_check (tau))
+      {
+        Error ("Radiation:sane_check CHECKING ff=%e, comp=%e, ind_comp=%e\n",
+  	     frac_ff, frac_comp, frac_ind_comp);
+      }
+  /* Calculate the reduction in the w of the photon bundle along with the average
+     weight in the cell */
+
+    if (tau > 0.0001)
+      {				/* Need differentiate between thick and thin cases */
+        x = exp (-tau);
+        p->w = w_out = w_in * x;
+        w_ave = (w_in - w_out) / tau;
+      }
+    else
+      {
+        tau2 = tau * tau;
+        p->w = w_out = w_in * (1. - tau + 0.5 * tau2);	/*Calculate to second order */
+        w_ave = w_in * (1. - 0.5 * tau + 0.1666667 * tau2);
+      }
+
+
+
+
+
+
+
+
   /*74a_ksl: 121215 -- Added to check on a problem photon */
   if (sane_check (p->w))
     {