Merge branch 'dev' of https://github.com/agnwinds/python into dev

source/photon_gen.c

@@ -847,32 +847,58 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
   double t, tref, teff (), tdisk ();
   double log_g, gref, geff (), gdisk ();
   double dr, r;
+  double logdr,logrmin,logrmax,logr;
   double f, ltot;
+  double rmax_temp;
   double q1;
-  int nrings;
+  int nrings,i,icheck;
   int spectype;
   double emit, emittance_bb (), emittance_continuum ();
 
   /* Calculate the reference temperature and luminosity of the disk */
   tref = tdisk (m, mdot, rmin);
+
+
   gref = gdisk (m, mdot, rmin);
+
 
   /* Now compute the apparent luminosity of the disk.  This is not actually used
      to determine how annulae are set up.  It is just used to populate geo.ltot.
      It can change if photons hitting the disk are allowed to raise the temperature
    */
 
-  ltot = 0;
-  dr = (rmax - rmin) / STEPS;
-  for (r = rmin; r < rmax; r += dr)
+  logrmax=log(rmax);
+  logrmin=log(rmin);
+  logdr=(logrmax-logrmin)/STEPS;
+
+  for (nrings = 0; nrings < NRINGS; nrings++) //Initialise the structure
   {
-    t = teff (tref, (r + 0.5 * dr) / rmin);
-    ltot += t * t * t * t * (2. * r + dr);
+    disk.nphot[nrings] = 0;
+    disk.nphot[nrings] = 0;
+	disk.r[nrings] = 0;
+	disk.t[nrings] = 0;	
+    disk.nhit[nrings] = 0;
+    disk.heat[nrings] = 0;
+    disk.ave_freq[nrings] = 0;
+    disk.w[nrings] = 0;
+    disk.t_hit[nrings] = 0;
   }
-  geo.lum_disk_init=ltot *= 2. * STEFAN_BOLTZMANN * PI * dr;
 
 
 
+
+  ltot = 0;
+
+ for (logr=logrmin;logr<logrmax;logr+=logdr)
+ {
+    r=exp(logr);
+    dr=exp(logr+logdr)-r;
+    t = teff (tref, (r + 0.5 * dr) / rmin);
+    ltot += t * t * t * t * (2. * r + dr) *dr;
+  }
+  geo.lum_disk_init=ltot *= 2. * STEFAN_BOLTZMANN * PI;
+
+
   /* Now establish the type of spectrum to create */
 
   if (ioniz_or_final == 1)
@@ -886,11 +912,16 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
    The extra factor of two arises because the disk radiates from both of its sides.
    */
 
-  q1 = 2. * PI * dr;
+  q1 = 2. * PI;
 
   (*ftot) = 0;
-  for (r = rmin; r < rmax; r += dr)
+  icheck=0;
+
+
+  for (logr=logrmin;logr<logrmax;logr+=logdr)
   {
+    r=exp(logr);
+    dr=exp(logr+logdr)-r;
     t = teff (tref, (r + 0.5 * dr) / rmin);
     log_g = log10 (geff (gref, (r + 0.5 * dr) / rmin));
 
@@ -901,12 +932,13 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
     else
     {
       emit = emittance_bb (freqmin, freqmax, t);
+
     }
-
-    (*ftot) += emit * (2. * r + dr);
+    (*ftot) += emit * (2. * r + dr) * dr;
   }
 
   (*ftot) *= q1;
+
 
 
   /* If *ftot is 0 in this energy range then all the photons come elsewhere, e. g. the star or BL  */
@@ -925,8 +957,12 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
   disk.v[0] = sqrt (G * geo.mstar / rmin);
   nrings = 1;
   f = 0;
-  for (r = rmin; r < rmax; r += dr)
+
+  i=0;
+  for (logr=logrmin;logr<logrmax;logr+=logdr)
   {
+    r=exp(logr);
+    dr=exp(logr+logdr)-r;
     t = teff (tref, (r + 0.5 * dr) / rmin);
     log_g = log10 (geff (gref, (r + 0.5 * dr) / rmin));
 
@@ -939,21 +975,21 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
       emit = emittance_bb (freqmin, freqmax, t);
     }
 
-    f += q1 * emit * (2. * r + dr);
+    f += q1 * emit * (2. * r + dr) * dr;
+	i++;
     /* EPSILON to assure that roundoffs don't affect result of if statement */
     if (f / (*ftot) * (NRINGS - 1) >= nrings)
     {
-      if (r <= disk.r[nrings - 1])
+      if (r <= disk.r[nrings - 1])  //If the radius we have reached is smaller than or equal to the last assigned radius - we make a tiny annulus
       {
         r = disk.r[nrings - 1] * (1. + 1.e-10);
       }
-
       disk.r[nrings] = r;
       disk.v[nrings] = sqrt (G * geo.mstar / r);
       nrings++;
       if (nrings >= NRINGS)
       {
-        Error_silent ("disk_init: Got to ftot %e at r %e < rmax %e. OK if freqs are high\n", f, r, rmax);
+//        Error_silent ("disk_init: Got to ftot %e at r %e < rmax %e. OK if freqs are high\n", f, r, rmax);		Not *really* an error, the error below deals with a *real* problem.		
         break;
       }
     }
@@ -964,16 +1000,16 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
     exit (0);
   }
 
-
-  disk.r[NRINGS - 1] = rmax;
-  disk.v[NRINGS - 1] = sqrt (G * geo.mstar / rmax);
+ 
+  disk.r[NRINGS - 1] = exp(logrmax); 
+  disk.v[NRINGS - 1] = sqrt (G * geo.mstar / disk.r[NRINGS - 1]);
 
 
   /* Now calculate the temperature and gravity of the annulae */
 
   for (nrings = 0; nrings < NRINGS - 1; nrings++)
   {
-    r = 0.5 * (disk.r[nrings + 1] + disk.r[nrings]);
+    r = 0.5 * (disk.r[nrings + 1] + disk.r[nrings]);	
     disk.t[nrings] = teff (tref, r / rmin);
     disk.g[nrings] = geff (gref, r / rmin);
   }
@@ -986,7 +1022,7 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_final, ftot)
     disk.heat[nrings] = 0;
     disk.ave_freq[nrings] = 0;
     disk.w[nrings] = 0;
-    disk.t_hit[nrings] = 0;
+    disk.t_hit[nrings] = 0;	
   }
 
   geo.lum_disk = ltot;
@@ -1152,6 +1188,8 @@ photo_gen_disk (p, weight, f1, f2, spectype, istart, nphot)
     p[i].freq /= (1. - dot (v, p[i].lmn) / C);
 
   }
+
+
   return (0);
 }
 

source/python.h

@@ -598,7 +598,7 @@ geo;
 
 
 /* xdisk is a structure that is used to store information about the disk in a system */
-#define NRINGS	301             /* The actual number of rings completely defined
+#define NRINGS	3001             /* The actual number of rings completely defined
                                    is NRINGS-1 ... or from 0 to NRINGS-2.  This is
                                    because you need an outer radius...but the rest
                                    of this element is not filled in. */

source/wind_updates2d.c

@@ -67,7 +67,7 @@ WindPtr (w);
 
   /*1108 NSH csum added to sum compton heating 1204 NSH icsum added to sum induced compton heating */
   double wtest, xsum, psum, fsum, lsum, csum, icsum, ausum;
-  double cool_sum,lum_sum; //1706 - the total cooling and luminosity of the wind
+  double cool_sum,lum_sum, rad_sum; //1706 - the total cooling and luminosity of the wind
   double apsum,aausum,abstot; //Absorbed photon energy from PI and auger
   double c_rec, n_rec, o_rec, fe_rec;   //1701- NSH more outputs to show cooling from a few other elements
   double c_lum, n_lum, o_lum, fe_lum;   //1708- NSH and luminosities as well
@@ -796,6 +796,13 @@ WindPtr (w);
   Log
     ("!!wind_update: Wind luminosity  %8.2e (recomb %8.2e ff %8.2e lines %8.2e) after update\n",
      lum_sum, geo.lum_rr, geo.lum_ff, geo.lum_lines);
+
+
+     rad_sum = wind_luminosity (xband.f1[0], xband.f2[xband.nbands-1]);       /*and we also call wind_luminosity to get the luminosities */
+
+    Log
+       ("!!wind_update: Rad  luminosity  %8.2e (recomb %8.2e ff %8.2e lines %8.2e) after update\n",
+         rad_sum, geo.lum_rr, geo.lum_ff, geo.lum_lines);
 
   Log
     ("!!wind_update: Wind cooling     %8.2e (recomb %8.2e ff %8.2e compton %8.2e DR %8.2e DI %8.2e lines %8.2e adiabatic %8.2e) after update\n",