text: #include <simtbx/nanoBragg/nanoBragg.h>
  
  
  //Contributed by James Holton, UCSF,LBNL,SLAC.
  
  namespace simtbx {
  namespace nanoBragg {
  
  
  /* constructor that takes a dxtbx "panel" detector model */
  nanoBragg::nanoBragg(
      const dxtbx::model::Detector& detector, // no default
      const dxtbx::model::Beam& beam,  // no default
      int verbose, // =0
      int panel_id) //=0
  {
      double temp;
      vec3 xyz;
  
      /* initialize to sensible default values */
      init_defaults();
  
      /* store arg verbosity as member */
      this->verbose=verbose;
  
      if(verbose) std::cout << detector << std::endl;
      if(verbose) std::cout << beam << std::endl;
  
      /* BEAM properties first */
  
      /* direction in 3-space of beam vector */
      xyz = beam.get_unit_s0();
      beam_vector[1] = xyz[0];
      beam_vector[2] = xyz[1];
      beam_vector[3] = xyz[2];
      unitize(beam_vector,beam_vector);
  
      /* central wavelength, in Angstrom */
      lambda0 = beam.get_wavelength()*1e-10;
  
      /* divergence, what are the DXTBX units? */
      temp = beam.get_divergence();
      if(temp>0.0) hdivrange = vdivrange = temp;
  
      /* assume this is photons/s, unless it is zero */
      temp = beam.get_flux();
      if(temp>0.0) flux = temp;
  
      /* assume this is Kahn polarization parameter */
      temp = beam.get_polarization_fraction();
      if(temp>=-1.0 && temp<=1.0) polarization = temp;
  
      /* dxtbx polarization points down B vector, we want the E vector */
      xyz = beam.get_polarization_normal();
      vert_vector[1] = xyz[0];
      vert_vector[2] = xyz[1];
      vert_vector[3] = xyz[2];
      unitize(vert_vector,vert_vector);
      cross_product(beam_vector,vert_vector,polar_vector);
      unitize(polar_vector,polar_vector);
  
  
      /* DETECTOR properties */
  
      /* size of the pixels in meters */
      pixel_size = detector[panel_id].get_pixel_size()[0]/1000.;
  
      /* pixel count in short and fast-axis directions */
      spixels = detector[panel_id].get_image_size()[1];
      fpixels = detector[panel_id].get_image_size()[0];
  
      /* allocate all the other arrays */
      init_detector();
  
      /* direction in 3-space of detector axes */
      beam_convention = CUSTOM;
      /* typically: 1 0 0 */
      fdet_vector[1] = detector[panel_id].get_fast_axis()[0];
      fdet_vector[2] = detector[panel_id].get_fast_axis()[1];
      fdet_vector[3] = detector[panel_id].get_fast_axis()[2];
      unitize(fdet_vector,fdet_vector);
      /* typically: 0 -1 0 */
      sdet_vector[1] = detector[panel_id].get_slow_axis()[0];
      sdet_vector[2] = detector[panel_id].get_slow_axis()[1];
      sdet_vector[3] = detector[panel_id].get_slow_axis()[2];
      unitize(sdet_vector,sdet_vector);
      /* set orthogonal vector to the detector pixel array */
      cross_product(fdet_vector,sdet_vector,odet_vector);
      unitize(odet_vector,odet_vector);
  
      /* dxtbx origin is location of outer corner of the first pixel */
      pix0_vector[1] = detector[panel_id].get_origin()[0]/1000.0;
      pix0_vector[2] = detector[panel_id].get_origin()[1]/1000.0;
      pix0_vector[3] = detector[panel_id].get_origin()[2]/1000.0;
      /* what is the point of closest approach between sample and detector? */
      Fclose = Xclose = -dot_product(pix0_vector,fdet_vector);
      Sclose = Yclose = -dot_product(pix0_vector,sdet_vector);
      close_distance = distance =  dot_product(pix0_vector,odet_vector);
  
      /* set beam centre */
      scitbx::vec2<double> dials_bc = detector[panel_id].get_beam_centre(beam.get_s0());
      Xbeam = dials_bc[0]/1000.0;
      Ybeam = dials_bc[1]/1000.0;
  
      /* detector sensor layer properties */
      detector_thick   = detector[panel_id].get_thickness();
      temp = detector[panel_id].get_mu();        // is this really a mu? or mu/rho ?
      if(temp>0.0) detector_attnlen = 1.0/temp;
  
      /* quantum_gain = amp_gain * electrooptical_gain, does not include capture_fraction */
      quantum_gain = detector[panel_id].get_gain();
  
      //adc_offset = detector[panel_id].ADC_OFFSET;
  
      /* SPINDLE properties */
  
      /* By default align the rotation axis with the detector fast direction */
      spindle_vector[1] = fdet_vector[1];
      spindle_vector[2] = fdet_vector[2];
      spindle_vector[3] = fdet_vector[3];
      unitize(spindle_vector,spindle_vector);
  
      /* NOT IMPLEMENTED: read in any other stuff?  */
      /*TODO: consider reading in a crystal model as well, showing params without crystal model can be confusing*/
      if (verbose)
        show_params();
  
      /* sensible initialization of all unititialized values */
      reconcile_parameters();
      SCITBX_ASSERT(distance > 0);  // safe to check,  especially with multi panel dxtbx models
  }
  
  
  // constructor for the nanoBragg class that takes most any member as an argument, defaults in nanoBragg_ext.cpp
  nanoBragg::nanoBragg(
          scitbx::vec2<int> detpixels_slowfast, // = 1024, 1024
          scitbx::vec3<int> Ncells_abc, // 1 1 1
          cctbx::uctbx::unit_cell unitcell, // lysozyme
          vec3 missets_deg, // 0 0 0
          vec2 beam_center_mm, // NAN NAN
          double distance_mm, // =100,
          double pixel_size_mm, // =0.1,
          double wavelength_A, // =1,
          double divergence_mrad, // =0,
          double dispersion_pct, // =0,
          double mosaic_spread_deg, // =-1
          int oversample, // =0 =auto
          int verbose)    // =0
  {
      /* initialize to sensible default values */
      init_defaults();
  
      /* now incorporate arguments to the constructor */
      this->spixels=detpixels_slowfast[0];
      this->fpixels=detpixels_slowfast[1];
      this->Na=Ncells_abc[0]; this->Nb=Ncells_abc[1]; this->Nc=Ncells_abc[2];
      this->a_A[0]=unitcell.parameters()[0];      // units = Angstrom
      this->b_A[0]=unitcell.parameters()[1];
      this->c_A[0]=unitcell.parameters()[2];
      this->alpha=unitcell.parameters()[3]/RTD;   // input units = Degrees
      this->beta=unitcell.parameters()[4]/RTD;
      this->gamma=unitcell.parameters()[5]/RTD;
      this->user_cell=true;
      this->misset[1]=missets_deg[0];
      this->misset[2]=missets_deg[1];
      this->misset[3]=missets_deg[2];
      if(this->misset[1] != 0.0 && this->misset[2] != 0.0 && this->misset[3] != 0.0) this->misset[0]=1;
      if(! isnan(beam_center_mm[0]) && ! isnan(beam_center_mm[1])) {
          this->user_beam = true;
          this->Xbeam = beam_center_mm[0]/1000.0;
          this->Ybeam = beam_center_mm[1]/1000.0;
      }
      this->distance=distance_mm/1000.0;
      this->pixel_size=pixel_size_mm/1000.0;
      this->lambda0=wavelength_A/1e10;
      this->hdivrange=divergence_mrad/1000.0;
      this->vdivrange=divergence_mrad/1000.0;
      this->dispersion=dispersion_pct/100.0;
      this->mosaic_spread=mosaic_spread_deg/RTD;
      this->oversample=oversample;
      if(oversample>0) this->user_oversample=true;
      this->verbose=verbose;
  
      /* sensible initialization of all unititialized values */
      reconcile_parameters();
  }
  
  
  /* sensible defaults */
  void
  nanoBragg::reconcile_parameters()
  {
      /* now that constructor arguments are applied, fill in blanks */
  
      /* allocate detector pixel array */
      init_detector();
  
      /* initialize fluence from flux and check sample size fits inside beam size */
      init_beam();
  
      /* initialize interpretation of detector and beam geometry */
      init_beamcenter();
  
      /* sanitize options for steps across divergence, dispersion, mosaic spread and spindle rotation */
      init_steps();
  
      /* reconcile different conventions of beam center input and detector position */
      update_beamcenter();
  
      /* set up interpolation if need be */
      init_interpolator();
  
      /* reconcile different types of input of cell or crystal orientation */
      init_cell();
  
      /* decide if oversampling is needed based on crystal volume and pixel size */
      update_oversample();
  
      /* read in structure factors */
      init_Fhkl();
  
      /* read in centrosymmetric background profile (not implemented) */
      init_background();
  
      /* display all actual phi values */
      show_phisteps();
  
      /* display all detector layers */
      show_detector_thicksteps();
  
      /* read in or generate x-ray source properties */
      init_sources();
  
      /* allocate and set up mosaic domains */
      init_mosaicity();
  
  
      if(verbose) printf("CONSTRUCT!!! %d   x-ray beam: %f %f %f\n",(int) raw_pixels.size(),this->beam_vector[1],this->beam_vector[2],this->beam_vector[3]);
  }
  
  
  /* sensible defaults */
  void
  nanoBragg::init_defaults()
  {
      double test = NAN;
      if(! isnan(test)) {
          if(verbose) printf("ERROR! isnan(NAN) != %g\n",test);
          exit(9);
      };
  
      /* optional file stuff, to be removed eventually? */
      matfilename = NULL;
      hklfilename = NULL;
      stolfilename = NULL;
      imginfilename = NULL;
      maskfilename = NULL;
  //    stoloutfilename = "output.stol\0";
      sourcefilename = NULL;
  //    floatfilename = "floatimage.bin\0";
  //    intfilename = "intimage.img\0";
  //    pgmfilename = "image.pgm\0";
  //    noisefilename = "noiseimage.img\0";
      infile = NULL;
  //    FILE *outfile;  // = NULL;
  //    FILE *stoloutfile;  // = NULL;
  
      /* default progress meter stuff */
      progress_pixel=progress_pixels=0;
      progress_meter = 1;
      babble = 1;
      printout = 0;
      printout_spixel=printout_fpixel=-1;
      verbose = 1;
  #ifdef NANOBRAGG_HAVE_CUDA
      device_Id = 0; // for running the CUDA code
  #endif
  
      /* default x-ray beam properties */
  //    double beam_vector[4] = {0,1,0,0};  this->beam_vector = beam_vector;
      beam_vector[0] = 0;
      beam_vector[1] = 1;
      beam_vector[2] = 0;
      beam_vector[3] = 0;
      coherent = 0;
      far_source = 1;
      round_div = 1;
      lambda=1e-10;
      lambda_of = NULL;
      mosaic_spread=-1.0;
      mosaic_umats = NULL;
      dispersion=0.0;dispstep=-1.0;
      lambda0 = 1.0e-10;
      hdiv=0;hdivstep=-1.0; hdivrange= -1.0;
      vdiv=0;vdivstep=-1.0; vdivrange= -1.0;
      source_path=0;source_distance = 10.0;
      hdiv_tic=vdiv_tic=disp_tic=mos_tic=0;
      divsteps=hdivsteps=vdivsteps=dispsteps=-1;
      mosaic_domains=-1;
      weight=0;
      source=sources=0;
      source_X=source_Y=source_Z=source_I=source_lambda=NULL;
      allocated_sources=0;
  
      /* Thomson cross section (m^2) */
  //    r_e_sqr = 7.94079248018965e-30;
      /* default incident x-ray fluence in photons/m^2 - exactly cancels Thomson cross section */
      fluence = 125932015286227086360700780544.0;
      flux=0.0;exposure=1.0;beamsize=1e-4;
  
      /* default sample size stuff */
      N=1;
      Na=Nb=Nc=1.0;
      xtalsize_max=xtalsize_a=xtalsize_b=xtalsize_c=0.0;
      reciprocal_pixel_size=0.0;
  
      xtal_shape = SQUARE;
      fudge=1;
      xtal_size_x   = 0;             /* m */
      xtal_size_y   = 0;             /* m */
      xtal_size_z   = 0;             /* m */
      xtal_density    = 1.0e6;         /* g/m^3 */
      xtal_molecular_weight = 18.0;    /* g/mol */
      xtal_volume=0.0;xtal_molecules = 0.0;
      /* scale factor = F^2*r_e_sqr*fluence*Avogadro*volume*density/molecular_weight
                             m^2     ph/m^2  /mol      m^3   g/m^3    g/mol   */
  
      /* amorphous material properties for background */
      amorphous_sample_x = 0.0;
      amorphous_sample_y = 0.0;
      amorphous_sample_z = 0.0;
      amorphous_volume = 0.0;
      amorphous_molecular_weight = 18.0;
      amorphous_density = 1.0e6;   // 1 g/cm^3
      amorphous_molecules = 0;
      /* water F = 2.57 in forward direction */
  
      /* default detector stuff */
      pixel_size = 0.1e-3;
      fpixels=spixels=0;
      pixels=0;
      distance = 100.0e-3;
      detsize_f = 102.4e-3;
      detsize_s = 102.4e-3;
      detector_attnlen= 234e-6; // default to Si at 1 A x-ray wavelength?
      detector_thick=0.0;
      detector_thickstep=0.0;
      detector_thicksteps=-1;
  //    double fdet_vector[4]  = {0,0,0,1};  this->fdet_vector = fdet_vector;
      fdet_vector[0] = 0;
      fdet_vector[1] = 0;
      fdet_vector[2] = 0;
      fdet_vector[3] = 1;
  //    double sdet_vector[4]  = {0,0,-1,0}; this->sdet_vector  = sdet_vector;
      sdet_vector[0] = 0;
      sdet_vector[1] = 0;
      sdet_vector[2] = -1;
      sdet_vector[3] = 0;
  //    double odet_vector[4]  = {0,1,0,0};  this->odet_vector  = odet_vector;
      odet_vector[0] = 0;
      odet_vector[1] = 1;
      odet_vector[2] = 0;
      odet_vector[3] = 0;
  //    double pix0_vector[4]  = {0,0,0,0};  this->pix0_vector  = pix0_vector;
      pix0_vector[0] = 0;
      pix0_vector[1] = 0;
      pix0_vector[2] = 0;
      pix0_vector[3] = 0;
      detector_rotx=0.0;detector_roty=0.0;detector_rotz=0.0;
  //    double twotheta_axis[4] = {0,0,1,0}; this->twotheta_axis = twotheta_axis;
      twotheta_axis[0] = 0;
      twotheta_axis[1] = 0;
      twotheta_axis[2] = 1;
      twotheta_axis[3] = 0;
      detector_pivot = BEAM;
      beam_convention = MOSFLM;
      detector_twotheta = 0.0;
      airpath=omega_pixel=omega_Rsqr_pixel=omega_sum = 0.0;
      curved_detector = 0;
      point_pixel= 0;
      Xbeam=NAN;Ybeam=NAN;
      Fbeam=NAN;Sbeam=NAN;
      Fdet=Sdet=Odet=Fdet0=Sdet0=0.0;
      Xclose=NAN;Yclose=NAN;close_distance=NAN;
      Fclose=NAN;Sclose=NAN;
      ORGX=NAN;ORGY=NAN;
      adc_offset = 40.0;
  
      /* use these to remember "user" inputs */
      user_beam = false;
      user_distance =false;
      user_mosdomains = false;
  
      /* scattering vectors */
  //    double incident[4] = {0,0,0,0};
  //    double diffracted[4],diffracted0[4] = {0,0,0,0};
  //    double scattering[4];
      stol=twotheta=theta=0.0;
  
      /* diffraction geometry stuff */
      costwotheta=sintwotheta=psi=0;
      xd=yd=zd=xd0=yd0=zd0=0.0;
  //    double Ewald[4],Ewald0[4],relp[4];
      dmin=0;
      integral_form = 0; // experimental: use integral form of spots to avoid need for oversampling
  
      /* polarization stuff */
  //    double polar_vector[4] = {0,0,0,1}; this->polar_vector = polar_vector;
      polar_vector[0] = 0;
      polar_vector[1] = 0;
      polar_vector[2] = 0;
      polar_vector[3] = 1;
  //  vert_vector[4];
      polar=1.0;
      polarization=0.0;
      nopolar = 0;
  
      /* sampling */
      steps=0;
      roi_xmin=-1;roi_xmax=-1;roi_ymin=-1;roi_ymax=-1;
      oversample = -1;
      recommended_oversample=subS=subF=0;
      user_oversample=false;
      subpixel_size=0.0;
  
      /* spindle */
      phi0=0.0;phistep=-1.0;osc=-1.0;
      phisteps=-1;
  //    double spindle_vector[4] = {0,0,0,1}; this->spindle_vector = spindle_vector;
      spindle_vector[0] = 0;
      spindle_vector[1] = 0;
      spindle_vector[2] = 0;
      spindle_vector[3] = 1;
  
      /* structure factor representation */
      Fhkl = NULL;
      F000 = 0.0;
      default_F = 0.0;
  
      /* background stuff */
      default_Fbg = 0.0;
      stol_of = NULL;
      Fbg_of  = NULL;
      allocated_stols = 0;
      nearest=0;
      stol_file_mult=1.0e10;
      ignore_values=0;
      Fmap_pixel=false;
  
      /* radial median filter stuff */
      pixels_in = NULL;
      bin_of = NULL;
      bin_start = NULL;
      invalid_pixel = NULL;
  //    double median,mad,deviate,sign;
  //    double sum_arej,avg_arej,sumd_arej,rms_arej,rmsd_arej;
      diffimage = NULL;
      stolimage = NULL;
      Fimage = NULL;
  
      pythony_indices.clear();
      pythony_amplitudes.clear();
      pythony_stolFbg.clear();
      pythony_source_XYZ.clear();
  
      /* intensity stats */
      max_I = 0.0;
      max_I_x = 0.0; max_I_y = 0.0;
      photon_scale = 0.0;
      intfile_scale = 0.0;
      pgm_scale = 0.0;
      sumn = 0;
      overloads = 0;
  
      /* image file data */
      floatimage = NULL;
      maskimage = NULL;
      intimage = NULL;
      pgmimage = NULL;
      allocated_pixels=0;
  
      /* random number seeds */
      seed = -time((time_t *)0);
  //    if(verbose) printf("random number seed = %u\n",seed);
      mosaic_seed = 12345678;
      calib_seed = 123456789;
  
      /* point-spread function parameters */
      psf_type = FIBER;
      psf_fwhm = 80e-6;
      psf_radius = 0;
      readout_noise = 3.0;
      flicker_noise = 0.0;
      calibration_noise = 0.03;
      spot_scale = 1.0;
      quantum_gain = 1.0;
  
      /* interpolation defaults: auto-detect */
      interpolate = 2;
      i1=0;i2=0;i3=0;
      sub_Fhkl=NULL;
  
      /* unit cell stuff */
      user_cell = false;
      user_matrix = false;
  //    double a[4] = {0,0,0,0}; this->a = a;
  //    double b[4] = {0,0,0,0}; this->b = b;
  //    double c[4] = {0,0,0,0}; this->c = c;
      a_A[0]=0;a_A[1]=0;a_A[2]=0;a_A[3]=0;
      b_A[0]=0;b_A[1]=0;b_A[2]=0;b_A[3]=0;
      c_A[0]=0;c_A[1]=0;c_A[2]=0;c_A[3]=0;
      a[0]=0;a[1]=0;a[2]=0;a[3]=0;
      b[0]=0;b[1]=0;b[2]=0;b[3]=0;
      c[0]=0;c[1]=0;c[2]=0;c[3]=0;
      alpha=0.0;beta=0.0;gamma=0.0;
  //    double misset[4] = {0,0,0,0}; this->misset = misset;
      misset[0] = 0;
      misset[1] = 0;
      misset[2] = 0;
      misset[3] = 0;
      user_umat = false;
      umat[0]=1;umat[1]=0;umat[2]=0;
      umat[3]=0;umat[4]=1;umat[5]=0;
      umat[6]=0;umat[7]=0;umat[8]=1;
  
      /* special options */
  //    calculate_noise = 1;
  //    write_pgm = 1;
  //    binary_spots = false;
  }
  // end of init_defaults
  
  
  
  /* initialize detector size from mm or pixel specs */
  void
  nanoBragg::init_detector()
      {
      /* fill in blanks */
      if(fpixels) {
          detsize_f = pixel_size*fpixels;
      }
      if(spixels) {
          detsize_s = pixel_size*spixels;
      }
      if(verbose){
          printf("fpixels= %d detsize_f/pixel_size=%g roundoff=%g\n",fpixels,detsize_f/pixel_size,ceil(detsize_f/pixel_size-0.5));
      }
  
      /* sanity check */
      if(allocated_pixels==0 && invalid_pixel!=NULL) {
          printf("ERROR: problem with memory allocation! \n");
          exit(9);
      }
  
      /* initialize Pythony pixel array */
      this->fpixels = static_cast<int>(ceil(detsize_f/pixel_size-0.5));
      this->spixels = static_cast<int>(ceil(detsize_s/pixel_size-0.5));
      pixels = this->fpixels*this->spixels;
  
      if(pixels!=allocated_pixels) {
          /* actually allocate memory */
          if(verbose>6) printf("(re)allocating %d %ld-byte doubles for raw_pixels array\n",pixels,sizeof(double));
          raw_pixels = af::flex_double(af::flex_grid<>(spixels,fpixels));
  
          if(invalid_pixel!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte bool invalid_pixel at %p\n",pixels,sizeof(bool),invalid_pixel);
              free(invalid_pixel);
          }
          if(bin_of!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte unsigned int bin_of at %p\n",pixels,sizeof(unsigned int),bin_of);
              free(bin_of);
          }
          if(diffimage!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte double diffimage at %p\n",2*pixels,sizeof(double),diffimage);
              free(diffimage);
          }
          if(stolimage!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte double stolimage at %p\n",pixels,sizeof(double),stolimage);
              free(stolimage);
          }
          if(Fimage!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte double Fimage at %p\n",pixels,sizeof(double),Fimage);
              free(Fimage);
          }
          if(intimage!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte unsigned short int intimage at %p\n",pixels,sizeof(unsigned short int),intimage);
              free(intimage);
          }
          if(pgmimage!=NULL) {
              if(verbose>6) printf("freeing %d %ld-byte unsigned char invalid_pixel at %p\n",pixels,sizeof(unsigned char),invalid_pixel);
              free(pgmimage);
          }
          if(verbose>6) printf("allocating %d pixels 7 times for image flags\n",pixels);
          invalid_pixel  =         (bool *) calloc(pixels+10,sizeof(bool));
          bin_of         = (unsigned int *) calloc(pixels+10,sizeof(unsigned int));
          diffimage      =       (double *) calloc(2*pixels+10,sizeof(double));
          stolimage      =       (double *) calloc(pixels+10,sizeof(double));
          Fimage         =       (double *) calloc(pixels+10,sizeof(double));
          intimage = (unsigned short int *) calloc(pixels+10,sizeof(unsigned short int));
          pgmimage =      (unsigned char *) calloc(pixels+10,sizeof(unsigned char));
      }
  
      allocated_pixels = pixels;
  }
  // end of init_detector
  
  
  /* initialize fluence from flux and check sample/beam size */
  void
  nanoBragg::init_beam()
  {
      /* get fluence from flux */
      if(flux != 0.0 && exposure > 0.0 && beamsize >= 0){
          fluence = flux*exposure/beamsize/beamsize;
      }
      if(beamsize >= 0){
          if(beamsize < xtal_size_y){
              if(verbose) printf("WARNING: clipping xtal (%lg m high) with beam (%lg m)\n",xtal_size_y,beamsize);
              xtal_size_y = beamsize;
          }
          if(beamsize < xtal_size_z){
              if(verbose) printf("WARNING: clipping xtal (%lg m wide) with beam (%lg m)\n",xtal_size_z,beamsize);
              xtal_size_z = beamsize;
          }
          if(beamsize < amorphous_sample_y){
              if(verbose) printf("WARNING: clipping amorphous sample (%lg m high) with beam (%lg m)\n",amorphous_sample_y,beamsize);
              amorphous_sample_y = beamsize;
          }
          if(beamsize < amorphous_sample_z){
              if(verbose) printf("WARNING: clipping amorphous sample (%lg m wide) with beam (%lg m)\n",amorphous_sample_z,beamsize);
              amorphous_sample_z = beamsize;
          }
      }
      if(exposure > 0.0)
      {
          flux = fluence/exposure*beamsize*beamsize;
      }
      /* straighten up sample properties */
  //    volume = sample_x*sample_y*sample_z;
  //    molecules = volume*density*Avogadro/molecular_weight;
  }
  // end of init_beam
  
  
  /* set up any uninitialized beam centers from others provided */
  void
  nanoBragg::init_beamcenter()
  {
      /* use XDS origin, if provided */
      if(! isnan(ORGX)) Xclose = ORGX*pixel_size;
      if(! isnan(ORGY)) Yclose = ORGY*pixel_size;
      /* default to center of detector */
      if(isnan(Xclose)) Xclose = (detsize_f - pixel_size)/2.0;
      if(isnan(Yclose)) Yclose = (detsize_s + pixel_size)/2.0;
      /* apply defaults for anything not initialized */
      if(isnan(Xbeam)) Xbeam = Xclose;
      if(isnan(Ybeam)) Ybeam = Yclose;
      if(isnan(Fbeam)) Fbeam = Fclose;
      if(isnan(Sbeam)) Sbeam = Sclose;
      if(roi_xmin < 0) roi_xmin = 0;
      if(roi_xmax < 0) roi_xmax = fpixels;
      if(roi_ymin < 0) roi_ymin = 0;
      if(roi_ymax < 0) roi_ymax = spixels;
      progress_pixels = (roi_xmax-roi_xmin+1)*(roi_ymax-roi_ymin+1);
  
      if(beam_convention == ADXV)
      {
          /* first pixel is at 0,0 pix and pixel_size,pixel_size*npixels mm */
          if(! user_beam)
          {
              Xbeam = (detsize_f + pixel_size)/2.0;
              Ybeam = (detsize_s - pixel_size)/2.0;
          }
             beam_vector[1]=  0;    beam_vector[2]=  0;    beam_vector[3]=  1;
             fdet_vector[1]=  1;    fdet_vector[2]=  0;    fdet_vector[3]=  0;
             sdet_vector[1]=  0;    sdet_vector[2]= -1;    sdet_vector[3]=  0;
             odet_vector[1]=  0;    odet_vector[2]=  0;    odet_vector[3]=  1;
           twotheta_axis[1]= -1;  twotheta_axis[2]=  0;  twotheta_axis[3]=  0;
            polar_vector[1]=  1;   polar_vector[2]=  0;   polar_vector[3]=  0;
          spindle_vector[1]=  1; spindle_vector[2]=  0; spindle_vector[3]=  0;
          Fbeam = Xbeam;
          Sbeam = detsize_s - Ybeam;
          detector_pivot = BEAM;
          if(verbose) printf("beam center in adxv convention\n");
      }
      if(beam_convention == MOSFLM)
      {
          /* first pixel is at 0.5,0.5 pix and pixel_size/2,pixel_size/2 mm */
          if(! user_beam)
          {
              Xbeam = (detsize_s + pixel_size)/2.0;
              Ybeam = (detsize_f + pixel_size)/2.0;
          }
             beam_vector[1]=  1;    beam_vector[2]=  0;    beam_vector[3]=  0;
             odet_vector[1]=  1;    odet_vector[2]=  0;    odet_vector[3]=  0;
             fdet_vector[1]=  0;    fdet_vector[2]=  0;    fdet_vector[3]=  1;
             sdet_vector[1]=  0;    sdet_vector[2]= -1;    sdet_vector[3]=  0;
           twotheta_axis[1]=  0;  twotheta_axis[2]=  0;  twotheta_axis[3]= -1;
            polar_vector[1]=  0;   polar_vector[2]=  0;   polar_vector[3]=  1;
          spindle_vector[1]=  0; spindle_vector[2]=  0; spindle_vector[3]=  1;
          Fbeam = Ybeam + 0.5*pixel_size;
          Sbeam = Xbeam + 0.5*pixel_size;
          detector_pivot = BEAM;
          if(verbose) printf("beam center in mosflm convention\n");
      }
      if(beam_convention == DENZO)
      {
          if(! user_beam)
          {
              Xbeam = (detsize_s + pixel_size)/2.0;
              Ybeam = (detsize_f + pixel_size)/2.0;
          }
             beam_vector[1]=  1;    beam_vector[2]=  0;    beam_vector[3]=  0;
             odet_vector[1]=  1;    odet_vector[2]=  0;    odet_vector[3]=  0;
             fdet_vector[1]=  0;    fdet_vector[2]=  0;    fdet_vector[3]=  1;
             sdet_vector[1]=  0;    sdet_vector[2]= -1;    sdet_vector[3]=  0;
           twotheta_axis[1]=  0;  twotheta_axis[2]=  0;  twotheta_axis[3]= -1;
            polar_vector[1]=  0;   polar_vector[2]=  0;   polar_vector[3]=  1;
          spindle_vector[1]=  0; spindle_vector[2]=  0; spindle_vector[3]=  1;
          Fbeam = Ybeam + 0.0*pixel_size;
          Sbeam = Xbeam + 0.0*pixel_size;
          detector_pivot = BEAM;
          if(verbose) printf("beam center in denzo convention\n");
      }
      if(beam_convention == XDS)
      {
          if(! user_beam)
          {
              Xbeam = Xclose;
              Ybeam = Yclose;
          }
             beam_vector[1]=  0;    beam_vector[2]=  0;    beam_vector[3]=  1;
             fdet_vector[1]=  1;    fdet_vector[2]=  0;    fdet_vector[3]=  0;
             sdet_vector[1]=  0;    sdet_vector[2]=  1;    sdet_vector[3]=  0;
             odet_vector[1]=  0;    odet_vector[2]=  0;    odet_vector[3]=  1;
           twotheta_axis[1]=  1;  twotheta_axis[2]=  0;  twotheta_axis[3]=  0;
            polar_vector[1]=  1;   polar_vector[2]=  0;   polar_vector[3]=  0;
          spindle_vector[1]=  1; spindle_vector[2]=  0; spindle_vector[3]=  0;
          Fbeam = Xbeam;
          Sbeam = Ybeam;
          detector_pivot = SAMPLE;
          if(verbose) printf("beam center in XDS convention\n");
      }
      if(beam_convention == DIALS)
      {
          if(! user_beam)
          {
              Xbeam = Xclose;
              Ybeam = Yclose;
          }
             beam_vector[1]=  0;    beam_vector[2]=  0;    beam_vector[3]= -1;
             fdet_vector[1]=  1;    fdet_vector[2]=  0;    fdet_vector[3]=  0;
             sdet_vector[1]=  0;    sdet_vector[2]= -1;    sdet_vector[3]=  0;
             odet_vector[1]=  0;    odet_vector[2]=  0;    odet_vector[3]= -1;
           twotheta_axis[1]=  1;  twotheta_axis[2]=  0;  twotheta_axis[3]=  0;
            polar_vector[1]=  1;   polar_vector[2]=  0;   polar_vector[3]=  0;
          spindle_vector[1]=  1; spindle_vector[2]=  0; spindle_vector[3]=  0;
          Fbeam = Xbeam;
          Sbeam = Ybeam;
          detector_pivot = SAMPLE;
          if(verbose) printf("beam center in DIALS convention\n");
      }
      if(beam_convention == CUSTOM)
      {
          if(! user_beam)
          {
              Xbeam = Xclose;
              Ybeam = Yclose;
          }
          Fbeam = Xbeam;
          Sbeam = Ybeam;
          Fclose = Xbeam;
          Sclose = Ybeam;
          if(verbose) printf("beam center in custom convention\n");
      }
  
      /* straighten up vectors */
      unitize(beam_vector,beam_vector);
      unitize(fdet_vector,fdet_vector);
      unitize(sdet_vector,sdet_vector);
      if(unitize(odet_vector,odet_vector) != 1.0)
      {
          if(verbose) printf("WARNING: auto-generating odet_vector\n");
          cross_product(fdet_vector,sdet_vector,odet_vector);
          unitize(odet_vector,odet_vector);
      }
      unitize(polar_vector,polar_vector);
      unitize(spindle_vector,spindle_vector);
      cross_product(beam_vector,polar_vector,vert_vector);
      unitize(vert_vector,vert_vector);
  
  }
  // init_beamcenter
  
  
  
  /* count up and sanitize steps across divergence, dispersion, mosaic spread and spindle rotation */
  void
  nanoBragg::init_steps()
  {
      /* default sampling logic */
      if(phisteps < 0){
          /* auto-select number of phi steps */
          if(osc < 0.0) {
              /* auto-select osc range */
              if(phistep <= 0.0) {
                  /* user doesn't care about anything */
                  phisteps = 1;
                  osc = 0.0;
                  phistep = 0.0;
              } else {
                  /* user doesn't care about osc or steps, but specified step */
                  osc = phistep;
                  phisteps = 2;
              }
          } else {
              /* user-speficied oscillation */
              if(phistep <= 0.0) {
                  /* osc specified, but nothing else */
                  phisteps = 2;
                  phistep = osc/2.0;
              } else {
                  /* osc and phi step specified */
                  phisteps = static_cast<int>(ceil(osc/phistep));
              }
          }
      } else {
          /* user-specified number of phi steps */
          if(phisteps == 0) phisteps = 1;
          if(osc < 0.0) {
              /* auto-select osc range */
              if(phistep <= 0.0) {
                  /* user cares only about number of steps */
                  osc = 1.0/RTD;
                  phistep = osc/phisteps;
              } else {
                  /* user doesn't care about osc, but specified step */
                  osc = phistep;
                  phisteps = 2;
              }
          } else {
              /* user-speficied oscillation */
              if(phistep < 0.0) {
                  /* osc and steps specified */
                  phistep = osc/phisteps;
              } else {
                  /* everything specified */
              }
          }
      }
  
      if(hdivsteps <= 0){
          /* auto-select number of steps */
          if(hdivrange < 0.0) {
              /* auto-select range */
              if(hdivstep <= 0.0) {
                  /* user doesn't care about anything */
                  hdivsteps = 1;
                  hdivrange = 0.0;
                  hdivstep = 0.0;
              } else {
                  /* user specified stepsize and nothing else */
                  hdivrange = hdivstep;
                  hdivsteps = 2;
              }
          } else {
              /* user-speficied range */
              if(hdivstep <= 0.0) {
                  /* range specified, but nothing else */
                  hdivstep = hdivrange;
                  hdivsteps = 2;
              } else {
                  /* range and step specified, but not number of steps */
                  hdivsteps = static_cast<int>(ceil(hdivrange/hdivstep));
              }
          }
      } else {
          /* user-specified number of steps */
          if(hdivrange < 0.0) {
              /* auto-select range */
              if(hdivstep <= 0.0) {
                  /* user cares only about number of steps */
                  hdivrange = 1.0;
                  hdivstep = hdivrange/hdivsteps;
              } else {
                  /* user doesn't care about range */
                  hdivrange = hdivstep;
                  hdivsteps = 2;
              }
          } else {
              /* user-speficied range */
              if(hdivstep <= 0.0) {
                  /* range and steps specified */
                  if(hdivsteps <=1 ) hdivsteps = 2;
                  hdivstep = hdivrange/(hdivsteps-1);
              } else {
                  /* everything specified */
              }
          }
      }
  
      if(vdivsteps <= 0){
          /* auto-select number of steps */
          if(vdivrange < 0.0) {
              /* auto-select range */
              if(vdivstep <= 0.0) {
                  /* user doesn't care about anything */
                  vdivsteps = 1;
                  vdivrange = 0.0;
                  vdivstep = 0.0;
              } else {
                  /* user specified stepsize and nothing else */
                  vdivrange = vdivstep;
                  vdivsteps = 2;
              }
          } else {
              /* user-speficied range */
              if(vdivstep <= 0.0) {
                  /* range specified, but nothing else */
                  vdivstep = vdivrange;
                  vdivsteps = 2;
              } else {
                  /* range and step specified, but not number of steps */
                  vdivsteps = static_cast<int>(ceil(vdivrange/vdivstep));
              }
          }
      } else {
          /* user-specified number of steps */
          if(vdivrange < 0.0) {
              /* auto-select range */
              if(vdivstep <= 0.0) {
                  /* user cares only about number of steps */
                  vdivrange = 1.0;
                  vdivstep = vdivrange/vdivsteps;
              } else {
                  /* user doesn't care about range */
                  vdivrange = vdivstep;
                  vdivsteps = 2;
              }
          } else {
              /* user-speficied range */
              if(vdivstep <= 0.0) {
                  /* range and steps specified */
                  if(vdivsteps <=1 ) vdivsteps = 2;
                  vdivstep = vdivrange/(vdivsteps-1);
              } else {
                  /* everything specified */
              }
          }
      }
  
  
      if(dispsteps <= 0){
          /* auto-select number of steps */
          if(dispersion < 0.0) {
              /* auto-select range */
              if(dispstep <= 0.0) {
                  /* user doesn't care about anything */
                  dispsteps = 1;
                  dispersion = 0.0;
                  dispstep = 0.0;
              } else {
                  /* user specified stepsize and nothing else */
                  dispersion = dispstep;
                  dispsteps = 2;
              }
          } else {
              /* user-speficied range */
              if(dispstep <= 0.0) {
                  /* range specified, but nothing else */
                  dispstep = dispersion;
                  dispsteps = 2;
              } else {
                  /* range and step specified, but not number of steps */
                  dispsteps = static_cast<int>(ceil(dispersion/dispstep));
              }
          }
      } else {
          /* user-specified number of steps */
          if(dispersion < 0.0) {
              /* auto-select range */
              if(dispstep <= 0.0) {
                  /* user cares only about number of steps */
                  dispersion = 1.0;
                  dispstep = dispersion/dispsteps;
              } else {
                  /* user doesn't care about range */
                  dispersion = dispstep;
                  dispsteps = 2;
              }
          } else {
              /* user-speficied range */
              if(dispstep <= 0.0) {
                  /* range and steps specified */
                  if(dispsteps <=1 ) dispsteps = 2;
                  dispstep = dispersion/(dispsteps-1);
              } else {
                  /* everything specified */
              }
          }
      }
  
  
      if(detector_thicksteps <= 0){
          /* auto-select number of steps */
          if(detector_thick < 0.0) {
              /* auto-select range */
              if(detector_thickstep <= 0.0) {
                  /* user doesn't care about anything */
                  detector_thicksteps = 1;
                  detector_thick = 0.0;
                  detector_thickstep = 0.0;
              } else {
                  /* user specified stepsize and nothing else */
                  detector_thick = detector_thickstep;
                  detector_thicksteps = 2;
              }
          } else {
              /* user-speficied range */
              if(detector_thickstep <= 0.0) {
                  /* range specified, but nothing else */
                  detector_thicksteps = 2;
                  detector_thickstep = detector_thick/detector_thicksteps;
              } else {
                  /* range and step specified, but not number of steps */
                  detector_thicksteps = ceil(detector_thick/detector_thickstep);
              }
          }
      } else {
          /* user-specified number of steps */
          if(detector_thick < 0.0) {
              /* auto-select range */
              if(detector_thickstep <= 0.0) {
                  /* user cares only about number of steps */
                  detector_thick = 0.5e-6;
                  detector_thickstep = detector_thick/detector_thicksteps;
              } else {
                  /* user doesn't care about range */
                  detector_thick = detector_thickstep;
                  detector_thicksteps = 2;
              }
          } else {
              /* user-speficied range */
              if(detector_thickstep <= 0.0) {
                  /* range and steps specified */
                  if(detector_thicksteps <=1 ) detector_thicksteps = 2;
                  detector_thickstep = detector_thick/(detector_thicksteps-1);
              } else {
                  /* everything specified */
              }
          }
      }
  
      if(mosaic_domains <= 0){
          /* auto-select number of domains */
          if(mosaic_spread < 0.0) {
              /* user doesn't care about anything */
              mosaic_domains = 1;
              mosaic_spread = 0.0;
          } else {
              /* user-speficied mosaicity, but not number of domains */
              if(mosaic_spread == 0.0)
              {
                  mosaic_domains = 1;
              }
              else
              {
                  if(verbose) printf("WARNING: finite mosaicity with only one domain! upping to 10 mosaic domains\n");
                  mosaic_domains = 10;
              }
          }
      } else {
          /* user-specified number of domains */
          if(mosaic_spread < 0.0) {
              /* number of domains specified, but no spread? */
              if(verbose) printf("WARNING: no mosaic spread specified.  setting mosaic_domains = 1\n");
              mosaic_spread = 0.0;
              mosaic_domains = 1;
          } else {
              /* user-speficied mosaicity and number of domains */
              if(mosaic_spread == 0.0)
              {
                  if(verbose) printf("WARNING: zero mosaic spread specified.  setting mosaic_domains = 1\n");
                  mosaic_domains = 1;
              }
          }
      }
  
  
      /* sanity checks */
      if(hdivrange <= 0.0 || hdivstep <= 0.0 || hdivsteps <= 0) {
          hdivsteps = 1;
          hdivrange = 0.0;
          hdivstep = 0.0;
      }
      if(vdivrange <= 0.0 || vdivstep <= 0.0 || vdivsteps <= 0) {
          vdivsteps = 1;
          vdivrange = 0.0;
          vdivstep = 0.0;
      }
      if(dispersion <= 0.0 || dispstep <= 0.0 || dispsteps <= 0) {
          dispsteps = 1;
          dispersion = 0.0;
          dispstep = 0.0;
      }
      if(detector_thick <= 0.0 || detector_thickstep <= 0.0 || detector_thicksteps <= 0) {
          detector_thicksteps = 1;
          detector_thick = 0.0;
          detector_thickstep = 0.0;
      }
      if(verbose) printf("dispersion= %g dispstep= %g  dispsteps= %d\n",dispersion,dispstep,dispsteps);
  
  }
  // end of init_steps()
  
  
  
  
  /* reconcile different conventions of beam center input and detector position */
  void
  nanoBragg::update_beamcenter()
  {
      /* initialize detector origin from a beam center and distance */
      /* there are different conventions here: mosflm, XDS etc */
  
      if(beam_convention == ADXV)
      {
          if(verbose) printf("adxv");
      }
      if(beam_convention == MOSFLM)
      {
          if(verbose) printf("mosflm");
      }
      if(beam_convention == XDS)
      {
          if(verbose) printf("xds");
      }
      if(beam_convention == DENZO)
      {
          if(verbose) printf("denzo");
      }
      if(beam_convention == DIALS)
      {
          if(verbose) printf("dials");
      }
      if(beam_convention == CUSTOM)
      {
          if(verbose) printf("custom");
      }
      if(verbose) printf(" convention selected.\n");
  
      /* first off, what is the relationship between the two "beam centers"? */
      rotate(odet_vector,vector,detector_rotx,detector_roty,detector_rotz);
      ratio = dot_product(beam_vector,vector);
      if(ratio == 0.0) {
          if(verbose) printf("WARNING: beam is parallel to detector surface! \n");
          ratio = DBL_MIN;
      }
      if(isnan(close_distance)) close_distance = fabs(ratio*distance);
      distance = close_distance/ratio;
  
      if(detector_pivot == SAMPLE){
          if(verbose) printf("pivoting detector around sample\n");
          /* initialize detector origin before rotating detector */
          pix0_vector[1] = -Fclose*fdet_vector[1]-Sclose*sdet_vector[1]+close_distance*odet_vector[1];
          pix0_vector[2] = -Fclose*fdet_vector[2]-Sclose*sdet_vector[2]+close_distance*odet_vector[2];
          pix0_vector[3] = -Fclose*fdet_vector[3]-Sclose*sdet_vector[3]+close_distance*odet_vector[3];
  
          /* now swing the detector origin around */
          rotate(pix0_vector,pix0_vector,detector_rotx,detector_roty,detector_rotz);
          rotate_axis(pix0_vector,pix0_vector,twotheta_axis,detector_twotheta);
      }
      /* now orient the detector plane */
      rotate(fdet_vector,fdet_vector,detector_rotx,detector_roty,detector_rotz);
      rotate(sdet_vector,sdet_vector,detector_rotx,detector_roty,detector_rotz);
      rotate(odet_vector,odet_vector,detector_rotx,detector_roty,detector_rotz);
  
      /* also apply orientation part of twotheta swing */
      rotate_axis(fdet_vector,fdet_vector,twotheta_axis,detector_twotheta);
      rotate_axis(sdet_vector,sdet_vector,twotheta_axis,detector_twotheta);
      rotate_axis(odet_vector,odet_vector,twotheta_axis,detector_twotheta);
  
      /* make sure beam center is preserved */
      if(detector_pivot == BEAM){
          if(verbose) printf("pivoting detector around direct beam spot\n");
          pix0_vector[1] = -Fbeam*fdet_vector[1]-Sbeam*sdet_vector[1]+distance*beam_vector[1];
          pix0_vector[2] = -Fbeam*fdet_vector[2]-Sbeam*sdet_vector[2]+distance*beam_vector[2];
          pix0_vector[3] = -Fbeam*fdet_vector[3]-Sbeam*sdet_vector[3]+distance*beam_vector[3];
      }
  
      /* what is the point of closest approach between sample and detector? */
      Fclose         = -dot_product(pix0_vector,fdet_vector);
      Sclose         = -dot_product(pix0_vector,sdet_vector);
      close_distance =  dot_product(pix0_vector,odet_vector);
  
      /* where is the direct beam now? */
      /* difference between beam impact vector and detector origin */
      newvector[1] = close_distance/ratio*beam_vector[1]-pix0_vector[1];
      newvector[2] = close_distance/ratio*beam_vector[2]-pix0_vector[2];
      newvector[3] = close_distance/ratio*beam_vector[3]-pix0_vector[3];
      /* extract components along detector vectors */
      Fbeam = dot_product(fdet_vector,newvector);
      Sbeam = dot_product(sdet_vector,newvector);
      distance = close_distance/ratio;
  
      /* find origin in XDS convention */
      ORGX=Fclose/pixel_size+0.5;
      ORGY=Sclose/pixel_size+0.5;
  
      /* find origin in DIALS convention
        from "CBF coordinate definition" at: https://sites.google.com/site/nexuscbf/home/cbf-dictionary
        The standard coordinate frame in imgCIF/CBF aligns the X-axis to the
        principal goniometer axis, chooses the Z-axis to point from the sample
        into the beam.  If the beam is not orthogonal to the X-axis, the Z-axis
        is the component of the vector points into the beam orthogonal to the
        X-axis.  The Y-axis is chosen to complete a right-handed axis system.
      */
      double cbfX[4],cbfY[4],cbfZ[4],negbeam[4];
      /* construct "X axis" */
      unitize(spindle_vector,cbfX);
      /* sample to source */
      vector_scale(beam_vector,negbeam,-1);
      /* construct "Z axis" */
      double comp = dot_product(negbeam,cbfX);
      cbfZ[1] = negbeam[1]-comp*cbfX[1];
      cbfZ[2] = negbeam[2]-comp*cbfX[2];
      cbfZ[3] = negbeam[3]-comp*cbfX[3];
      unitize(cbfZ,cbfZ);
      /* construct "Y axis" */
      cross_product(cbfZ,cbfX,cbfY);
      unitize(cbfY,cbfY);
      /* now get the components */
      dials_origin[1] = 1000.0*dot_product(pix0_vector,cbfX);
      dials_origin[2] = 1000.0*dot_product(pix0_vector,cbfY);
      dials_origin[3] = 1000.0*dot_product(pix0_vector,cbfZ);
  
      /* find the beam in the detector frame, for XDS */
      newvector[1] = dot_product(beam_vector,fdet_vector);
      newvector[2] = dot_product(beam_vector,sdet_vector);
      newvector[3] = dot_product(beam_vector,odet_vector);
      if(verbose) printf("XDS incident beam: %g %g %g\n",newvector[1],newvector[2],newvector[3]);
  }
  // end of update_beamcenter()
  
  
  
  
  /* automatically decide if we are interpolating, allocate memory if yes */
  void
  nanoBragg::init_interpolator()
  {
      /* free any previous allocations */
      if(sub_Fhkl != NULL) {
          for (h0=0; h0<=5;h0++) {
              for (k0=0; k0<=5;k0++) {
                  if(verbose>6) printf("freeing %d %ld-byte double Fhkl[%d][%d] at %p\n",5,sizeof(double),h0,k0,sub_Fhkl[h0][k0]);
                  free(*(*(sub_Fhkl +h0)+k0));
              }
              if(verbose>6) printf("freeing %d %ld-byte double* Fhkl[%d] at %p\n",5,sizeof(double*),h0,sub_Fhkl[h0]);
              free(*(sub_Fhkl +h0));
          }
          if(verbose>6) printf("freeing %d %ld-byte double** Fhkl at %p\n",5,sizeof(double**),sub_Fhkl);
          free(sub_Fhkl);
      }
  
      if(interpolate > 1){
          /* no user options */
          if(( Na <= 2) || (Nb <= 2) || (Nc <= 2)){
              if(verbose) printf("auto-selected tricubic interpolation of structure factors\n");
              interpolate = 1;
          }
          else
          {
              if(verbose) printf("auto-selected no interpolation\n");
              interpolate = 0;
          }
      }
      if(interpolate){
          /* allocate interpolation array */
          sub_Fhkl = (double***) calloc(6,sizeof(double**));
          for (h0=0; h0<=5;h0++) {
              *(sub_Fhkl +h0) = (double**) calloc(6,sizeof(double*));
              for (k0=0; k0<=5;k0++) {
                  *(*(sub_Fhkl +h0)+k0) = (double*) calloc(6,sizeof(double));
              }
          }
      }
  }
  // end of init_interpolator()
  
  
  
  
  /* compute A matrix from cell, missets, mosflm matrix file, etc. */
  void
  nanoBragg::init_cell()
  {
      /* do not run if things are not going to work */
      if(! user_cell && ! user_matrix && matfilename == NULL)
      {
          if(verbose) printf("ERROR: cannot initialize without a cell\n");
          return;
      }
      if(verbose>1) printf("raw CELL %f %f %f %f %f %f   %d\n",a_A[0],b_A[0],c_A[0],alpha,beta,gamma,user_cell);
      /* user-specified unit cell */
      if(user_cell)
      {
          /* a few random defaults */
          if(b_A[0]  <= 0.0) b_A[0] = a_A[0];
          if(c_A[0]  <= 0.0) c_A[0] = a_A[0];
          if(alpha <= 0.0) alpha = M_PI/2;
          if(beta  <= 0.0) beta  = M_PI/2;
          if(gamma <= 0.0) gamma = M_PI/2;
  
          /* get cell volume from angles */
          aavg = (alpha+beta+gamma)/2;
          skew = sin(aavg)*sin(aavg-alpha)*sin(aavg-beta)*sin(aavg-gamma);
          if(skew<0.0) skew=-skew;
          V_cell = 2.0*a_A[0]*b_A[0]*c_A[0]*sqrt(skew);
          if(V_cell <= 0.0)
          {
              if(verbose) printf("WARNING: impossible unit cell volume: %g\n",V_cell);
              V_cell = DBL_MIN;
          }
          V_star = 1.0/V_cell;
  
          /* now get reciprocal-cell lengths from the angles and volume */
          a_star[0] = b_A[0]*c_A[0]*sin(alpha)*V_star;
          b_star[0] = c_A[0]*a_A[0]*sin(beta)*V_star;
          c_star[0] = a_A[0]*b_A[0]*sin(gamma)*V_star;
          if(a_star[0] <= 0.0 || b_star[0] <= 0.0 || c_star[0] <= 0.0)
          {
              if(verbose) printf("WARNING: impossible reciprocal cell lengths: %g %g %g\n",
                  a_star[0],b_star[0],c_star[0]);
              /* make up something non-zero */
              a_star[0] = fabs(a_star[0]);
              b_star[0] = fabs(b_star[0]);
              c_star[0] = fabs(c_star[0]);
              if(a_star[0] <= 0.0) a_star[0] = DBL_MIN;
              if(b_star[0] <= 0.0) b_star[0] = DBL_MIN;
              if(c_star[0] <= 0.0) c_star[0] = DBL_MIN;
          }
  
          /* for fun, compute the reciprocal-cell angles from direct-cell angles */
          sin_alpha_star = a_A[0]*V_star/b_star[0]/c_star[0];
          sin_beta_star  = b_A[0]*V_star/a_star[0]/c_star[0];
          sin_gamma_star = c_A[0]*V_star/a_star[0]/b_star[0];
          cos_alpha_star = (cos(beta)*cos(gamma)-cos(alpha))/(sin(beta)*sin(gamma));
          cos_beta_star  = (cos(gamma)*cos(alpha)-cos(beta))/(sin(gamma)*sin(alpha));
          cos_gamma_star = (cos(alpha)*cos(beta)-cos(gamma))/(sin(alpha)*sin(beta));
          if(sin_alpha_star>1.0000001 || sin_alpha_star<-1.0000001 ||
             sin_beta_star >1.0000001 || sin_beta_star <-1.0000001 ||
             sin_gamma_star>1.0000001 || sin_gamma_star<-1.0000001 ||
             cos_alpha_star>1.0000001 || cos_alpha_star<-1.0000001 ||
             cos_beta_star >1.0000001 || cos_beta_star <-1.0000001 ||
             cos_gamma_star>1.0000001 || cos_gamma_star<-1.0000001 )
          {
              if(verbose) printf("WARNING: oddball reciprocal cell angles:\n");
              if(verbose) printf("sin(alpha_star) = %.25g\n",sin_alpha_star);
              if(verbose) printf("cos(alpha_star) = %.25g\n",cos_alpha_star);
              if(verbose) printf("sin(beta_star)  = %.25g\n",sin_beta_star);
              if(verbose) printf("cos(beta_star)  = %.25g\n",cos_beta_star);
              if(verbose) printf("sin(gamma_star) = %.25g\n",sin_gamma_star);
              if(verbose) printf("cos(gamma_star) = %.25g\n",cos_gamma_star);
          }
          if(sin_alpha_star>1.0) sin_alpha_star=1.0;
          if(sin_beta_star >1.0) sin_beta_star =1.0;
          if(sin_gamma_star>1.0) sin_gamma_star=1.0;
          if(sin_alpha_star<-1.0) sin_alpha_star=-1.0;
          if(sin_beta_star <-1.0) sin_beta_star =-1.0;
          if(sin_gamma_star<-1.0) sin_gamma_star=-1.0;
          if(cos_alpha_star*cos_alpha_star>1.0) cos_alpha_star=1.0;
          if(cos_beta_star *cos_beta_star >1.0) cos_beta_star=1.0;
          if(cos_gamma_star*cos_gamma_star>1.0) cos_gamma_star=1.0;
          alpha_star = atan2(sin_alpha_star,cos_alpha_star);
          beta_star  = atan2(sin_beta_star ,cos_beta_star );
          gamma_star = atan2(sin_gamma_star,cos_gamma_star);
  
          if(verbose>1) printf("CELL %f %f %f %f %f %f   %d\n",a_A[0],b_A[0],c_A[0],alpha,beta,gamma,user_cell);
  
          /* construct default orientation */
          a_star[1] = a_star[0];
          b_star[1] = b_star[0]*cos_gamma_star;
          c_star[1] = c_star[0]*cos_beta_star;
          a_star[2] = 0.0;
          b_star[2] = b_star[0]*sin_gamma_star;
          c_star[2] = c_star[0]*(cos_alpha_star-cos_beta_star*cos_gamma_star)/sin_gamma_star;
          a_star[3] = 0.0;
          b_star[3] = 0.0;
          c_star[3] = c_star[0]*V_cell/(a_A[0]*b_A[0]*c_A[0]*sin_gamma_star);
      }
  
      /* load the lattice orientation (reciprocal cell vectors) from a mosflm matrix */
      if(matfilename != NULL)
      {
          infile = fopen(matfilename,"r");
          if(infile != NULL)
          {
              if(verbose) printf("reading %s\n",matfilename);
              if(! fscanf(infile,"%lg%lg%lg",a_star+1,b_star+1,c_star+1)) {perror("fscanf");};
              if(! fscanf(infile,"%lg%lg%lg",a_star+2,b_star+2,c_star+2)) {perror("fscanf");};
              if(! fscanf(infile,"%lg%lg%lg",a_star+3,b_star+3,c_star+3)) {perror("fscanf");};
              fclose(infile);
  
              /* mosflm A matrix includes the wavelength, so remove it */
              /* calculate reciprocal cell lengths, store in 0th element */
              vector_scale(a_star,a_star,1e-10/lambda0);
              vector_scale(b_star,b_star,1e-10/lambda0);
              vector_scale(c_star,c_star,1e-10/lambda0);
          }
      }
  
      if(verbose>1) printf("MISSET %f %f %f \n",misset[1],misset[2],misset[3]);
      /* apply any user-provided unitary rotation matrix? */
      if(user_umat)
      {
          rotate_umat(a_star,a_star,umat);
          rotate_umat(b_star,b_star,umat);
          rotate_umat(c_star,c_star,umat);
      }
  
      if(verbose>1) printf("MISSET %f %f %f \n",misset[1],misset[2],misset[3]);
      /* apply any missetting angle */
      if(misset[0] > 0.0)
      {
          rotate(a_star,a_star,misset[1],misset[2],misset[3]);
          rotate(b_star,b_star,misset[1],misset[2],misset[3]);
          rotate(c_star,c_star,misset[1],misset[2],misset[3]);
      }
  
      /* various cross products */
      cross_product(a_star,b_star,a_star_cross_b_star);
      cross_product(b_star,c_star,b_star_cross_c_star);
      cross_product(c_star,a_star,c_star_cross_a_star);
  
      /* reciprocal lattice vector "a_star" is defined as perpendicular to both b and c, and must also preserve volume
         converse is true for direct-space lattice: a is perpendicular to both b_star and c_star
         a = ( b_star cross c_star ) / V_star    */
  
      /* reciprocal unit cell volume, but is it lambda-corrected? */
      V_star = dot_product(a_star,b_star_cross_c_star);
  
      /* make sure any user-supplied cell takes */
      if(user_cell)
      {
          /* a,b,c and V_cell were generated above */
  
          /* force the cross-product vectors to have proper magnitude: b_star X c_star = a*V_star */
          vector_rescale(b_star_cross_c_star,b_star_cross_c_star,a_A[0]/V_cell);
          vector_rescale(c_star_cross_a_star,c_star_cross_a_star,b_A[0]/V_cell);
          vector_rescale(a_star_cross_b_star,a_star_cross_b_star,c_A[0]/V_cell);
          V_star = 1.0/V_cell;
      }
  
      /* direct-space cell volume in A^3 */
      V_cell = 1.0/V_star;
  
      /* generate direct-space cell vectors, also updates magnitudes */
      vector_scale(b_star_cross_c_star,a_A,V_cell);
      vector_scale(c_star_cross_a_star,b_A,V_cell);
      vector_scale(a_star_cross_b_star,c_A,V_cell);
  
      /* now that we have direct-space vectors, re-generate the reciprocal ones */
      cross_product(a_A,b_A,a_cross_b);
      cross_product(b_A,c_A,b_cross_c);
      cross_product(c_A,a_A,c_cross_a);
      vector_scale(b_cross_c,a_star,V_star);
      vector_scale(c_cross_a,b_star,V_star);
      vector_scale(a_cross_b,c_star,V_star);
  
      /* for fun, calculate the cell angles too */
      sin_alpha = a_star[0]*V_cell/b_A[0]/c_A[0];
      sin_beta  = b_star[0]*V_cell/a_A[0]/c_A[0];
      sin_gamma = c_star[0]*V_cell/a_A[0]/b_A[0];
      cos_alpha = dot_product(b_A,c_A)/b_A[0]/c_A[0];
      cos_beta  = dot_product(a_A,c_A)/a_A[0]/c_A[0];
      cos_gamma = dot_product(a_A,b_A)/a_A[0]/b_A[0];
      if(sin_alpha>1.0000001 || sin_alpha<-1.0000001 ||
         sin_beta >1.0000001 || sin_beta <-1.0000001 ||
         sin_gamma>1.0000001 || sin_gamma<-1.0000001 ||
         cos_alpha>1.0000001 || cos_alpha<-1.0000001 ||
         cos_beta >1.0000001 || cos_beta <-1.0000001 ||
         cos_gamma>1.0000001 || cos_gamma<-1.0000001 )
      {
          if(verbose) printf("WARNING: oddball cell angles:\n");
              if(verbose) printf("sin_alpha = %.25g\n",sin_alpha);
              if(verbose) printf("cos_alpha = %.25g\n",cos_alpha);
              if(verbose) printf("sin_beta  = %.25g\n",sin_beta);
              if(verbose) printf("cos_beta  = %.25g\n",cos_beta);
              if(verbose) printf("sin_gamma = %.25g\n",sin_gamma);
              if(verbose) printf("cos_gamma = %.25g\n",cos_gamma);
      }
      if(sin_alpha>1.0) sin_alpha=1.0;
      if(sin_beta >1.0) sin_beta =1.0;
      if(sin_gamma>1.0) sin_gamma=1.0;
      if(sin_alpha<-1.0) sin_alpha=-1.0;
      if(sin_beta <-1.0) sin_beta =-1.0;
      if(sin_gamma<-1.0) sin_gamma=-1.0;
      if(cos_alpha*cos_alpha>1.0) cos_alpha=1.0;
      if(cos_beta *cos_beta >1.0) cos_beta=1.0;
      if(cos_gamma*cos_gamma>1.0) cos_gamma=1.0;
      alpha = atan2(sin_alpha,cos_alpha);
      beta  = atan2(sin_beta ,cos_beta );
      gamma = atan2(sin_gamma,cos_gamma);
  
  
      /* reciprocal cell angles */
      sin_alpha_star = a_A[0]*V_star/b_star[0]/c_star[0];
      sin_beta_star  = b_A[0]*V_star/a_star[0]/c_star[0];
      sin_gamma_star = c_A[0]*V_star/a_star[0]/b_star[0];
      cos_alpha_star = dot_product(b_star,c_star)/b_star[0]/c_star[0];
      cos_beta_star  = dot_product(a_star,c_star)/a_star[0]/c_star[0];
      cos_gamma_star = dot_product(a_star,b_star)/a_star[0]/b_star[0];
      if(sin_alpha_star>1.0000001 || sin_alpha_star<-1.0000001 ||
         sin_beta_star >1.0000001 || sin_beta_star <-1.0000001 ||
         sin_gamma_star>1.0000001 || sin_gamma_star<-1.0000001 ||
         cos_alpha_star>1.0000001 || cos_alpha_star<-1.0000001 ||
         cos_beta_star >1.0000001 || cos_beta_star <-1.0000001 ||
         cos_gamma_star>1.0000001 || cos_gamma_star<-1.0000001 )
      {
              if(verbose) printf("WARNING: oddball reciprocal cell angles:\n");
              if(verbose) printf("sin(alpha_star) = %.25g\n",sin_alpha_star);
              if(verbose) printf("cos(alpha_star) = %.25g\n",cos_alpha_star);
              if(verbose) printf("sin(beta_star)  = %.25g\n",sin_beta_star);
              if(verbose) printf("cos(beta_star)  = %.25g\n",cos_beta_star);
              if(verbose) printf("sin(gamma_star) = %.25g\n",sin_gamma_star);
              if(verbose) printf("cos(gamma_star) = %.25g\n",cos_gamma_star);
      }
      if(sin_alpha_star>1.0) sin_alpha_star=1.0;
      if(sin_beta_star >1.0) sin_beta_star =1.0;
      if(sin_gamma_star>1.0) sin_gamma_star=1.0;
      if(sin_alpha_star<-1.0) sin_alpha_star=-1.0;
      if(sin_beta_star <-1.0) sin_beta_star =-1.0;
      if(sin_gamma_star<-1.0) sin_gamma_star=-1.0;
      if(cos_alpha_star*cos_alpha_star>1.0) cos_alpha_star=1.0;
      if(cos_beta_star *cos_beta_star >1.0) cos_beta_star=1.0;
      if(cos_gamma_star*cos_gamma_star>1.0) cos_gamma_star=1.0;
      alpha_star = atan2(sin_alpha_star,cos_alpha_star);
      beta_star  = atan2(sin_beta_star ,cos_beta_star );
      gamma_star = atan2(sin_gamma_star,cos_gamma_star);
  
      if(verbose) printf("Unit Cell: %g %g %g %g %g %g\n", a_A[0],b_A[0],c_A[0],alpha*RTD,beta*RTD,gamma*RTD);
      if(verbose) printf("Recp Cell: %g %g %g %g %g %g\n", a_star[0],b_star[0],c_star[0],alpha_star*RTD,beta_star*RTD,gamma_star*RTD);
      if(verbose) printf("volume = %g A^3\n",V_cell);
  
      /* print out the real-space matrix */
      if(verbose) printf("real-space cell vectors (Angstrom):\n");
      if(verbose) printf("     %-10s  %-10s  %-10s\n","a","b","c");
      if(verbose) printf("X: %11.8f %11.8f %11.8f\n",a_A[1],b_A[1],c_A[1]);
      if(verbose) printf("Y: %11.8f %11.8f %11.8f\n",a_A[2],b_A[2],c_A[2]);
      if(verbose) printf("Z: %11.8f %11.8f %11.8f\n",a_A[3],b_A[3],c_A[3]);
      if(verbose) printf("reciprocal-space cell vectors (Angstrom^-1):\n");
      if(verbose) printf("     %-10s  %-10s  %-10s\n","a_star","b_star","c_star");
      if(verbose) printf("X: %11.8f %11.8f %11.8f\n",a_star[1],b_star[1],c_star[1]);
      if(verbose) printf("Y: %11.8f %11.8f %11.8f\n",a_star[2],b_star[2],c_star[2]);
      if(verbose) printf("Z: %11.8f %11.8f %11.8f\n",a_star[3],b_star[3],c_star[3]);
  
      /* now convert Angstrom to meters */
      vector_scale(a_A,a,1e-10);
      vector_scale(b_A,b,1e-10);
      vector_scale(c_A,c,1e-10);
  
      /* define phi=0 mosaic=0 crystal orientation */
      vector_scale(a,a0,1.0);
      vector_scale(b,b0,1.0);
      vector_scale(c,c0,1.0);
  
      /* define phi=0 crystal orientation */
      vector_scale(a,ap,1.0);
      vector_scale(b,bp,1.0);
      vector_scale(c,cp,1.0);
  }
  // end of init_cell()
  
  
  /* compute sample size and required oversampling */
  void
  nanoBragg::update_oversample()
  {
      /* now we know the cell, calculate crystal size in meters */
      if(xtal_size_x > 0) Na = ceil(xtal_size_x/a[0]-1e-6);
      if(xtal_size_y > 0) Nb = ceil(xtal_size_y/b[0]-1e-6);
      if(xtal_size_z > 0) Nc = ceil(xtal_size_z/c[0]-1e-6);
      if(Na <= 1.0) Na = 1.0;
      if(Nb <= 1.0) Nb = 1.0;
      if(Nc <= 1.0) Nc = 1.0;
      xtalsize_a = a[0]*Na;
      xtalsize_b = b[0]*Nb;
      xtalsize_c = c[0]*Nc;
      if(verbose) printf("crystal is %g x %g x %g microns\n",xtalsize_a*1e6,xtalsize_b*1e6,xtalsize_c*1e6);
      xtalsize_max = xtalsize_a;
      if(xtalsize_max < xtalsize_b) xtalsize_max = xtalsize_b;
      if(xtalsize_max < xtalsize_c) xtalsize_max = xtalsize_c;
      reciprocal_pixel_size = lambda0*distance/pixel_size;
      recommended_oversample = static_cast<int>(ceil(3.0 * xtalsize_max/reciprocal_pixel_size));
      if(recommended_oversample <= 0) recommended_oversample = 1;
      if(! user_oversample) {
          oversample = recommended_oversample;
          if(verbose) printf("auto-selected %d-fold oversampling\n",oversample);
      }
      if(oversample < recommended_oversample)
      {
          if(verbose)
          {
          printf("WARNING: maximum dimension of xtal is %g A\n",xtalsize_max*1e10);
          printf("         but reciprocal pixel size is %g A\n", reciprocal_pixel_size*1e10 );
          printf("         intensity may vary significantly across a pixel!\n");
              printf("         recommend oversample=%d to work around this\n",recommended_oversample);
          }
      }
  
      /* rough estimate of sample properties */
      xtal_size_x = xtalsize_a;
      xtal_size_y = xtalsize_b;
      xtal_size_z = xtalsize_c;
      xtal_volume = xtal_size_x*xtal_size_y*xtal_size_z;
      xtal_density = 1.2e6;
      xtal_molecules = Na*Nb*Nc;
      xtal_molecular_weight = xtal_volume*xtal_density*Avogadro/xtal_molecules;
      if(verbose) printf("approximate MW = %g\n",xtal_molecular_weight);
  
  }
  // end of update_oversample()
  
  
  
  
  /* read in structure factors vs hkl */
  void
  nanoBragg::init_Fhkl()
  {
      /* free any previous allocations */
      if(Fhkl != NULL) {
          for (h0=0; h0<=h_range;h0++) {
              for (k0=0; k0<=k_range;k0++) {
                  if(verbose>6) printf("freeing %d %ld-byte double Fhkl[%d][%d] at %p\n",l_range+1,sizeof(double),h0,k0,Fhkl[h0][k0]);
                  free(Fhkl[h0][k0]);
              }
              if(verbose>6) printf("freeing %d %ld-byte double* Fhkl[%d] at %p\n",k_range+1,sizeof(double*),h0,Fhkl[h0]);
              free(Fhkl[h0]);
          }
          if(verbose>6) printf("freeing %d %ld-byte double** Fhkl at %p\n",h_range+1,sizeof(double**),Fhkl);
          free(Fhkl);
      }
      hkls = 0;
  
      /* load the structure factors */
      if(hklfilename == NULL)
      {
          /* try to recover Fs from a previous run? */
          h_min = h_max = 0;
          k_min = k_max = 0;
          l_min = l_max = 0;
      }
      else
      {
          infile = fopen(hklfilename,"r");
          if(infile == NULL)
          {
              if(verbose) printf("ERROR: unable to open %s.",hklfilename);
              exit(9);
          }
          h_min=k_min=l_min=static_cast<int>(1e9);
          h_max=k_max=l_max=static_cast<int>(-1e9);
          if(verbose) printf("counting entries in %s\n",hklfilename);
          while(4 == fscanf(infile,"%lg%lg%lg%lg",&h,&k,&l,&F_cell)){
              if(verbose && h != ceil(h-0.4)) printf("WARNING: non-integer value for h (%g) at line %d\n",h,hkls);
              if(verbose && k != ceil(k-0.4)) printf("WARNING: non-integer value for k (%g) at line %d\n",k,hkls);
              if(verbose && l != ceil(l-0.4)) printf("WARNING: non-integer value for l (%g) at line %d\n",l,hkls);
              if(h_min > h) h_min = static_cast<int>(h);
              if(k_min > k) k_min = static_cast<int>(k);
              if(l_min > l) l_min = static_cast<int>(l);
              if(h_max < h) h_max = static_cast<int>(h);
              if(k_max < k) k_max = static_cast<int>(k);
              if(l_max < l) l_max = static_cast<int>(l);
              ++hkls;
          }
          rewind(infile);
      }
  
      if(pythony_indices.size())
      {
          if(verbose) printf(" noticed pythony_indices.size() = %ld\n",pythony_indices.size());
          /* need to know how much memory to allocate for Fhkl array */
          h_min=k_min=l_min=1e9;
          h_max=k_max=l_max=-1e9;
          miller_t hkl;
          for (i=0; i < pythony_indices.size(); ++i)
          {
              hkl = pythony_indices[i];
              if(pythony_amplitudes.size()) F_cell = pythony_amplitudes[i];
              h0 = hkl[0];
              k0 = hkl[1];
              l0 = hkl[2];
              if(verbose>9) printf("DEBUG: hkl# %d : %d %d %d = %g\n",i,h0,k0,l0,F_cell);
              if(h_min > h0) h_min = h0;
              if(k_min > k0) k_min = k0;
              if(l_min > l0) l_min = l0;
              if(h_max < h0) h_max = h0;
              if(k_max < k0) k_max = k0;
              if(l_max < l0) l_max = l0;
              ++hkls;
          }
      }
  
      if(1)
      {
          h_range = h_max - h_min + 1;
          k_range = k_max - k_min + 1;
          l_range = l_max - l_min + 1;
  
          if(verbose) printf("h: %d - %d\n",h_min,h_max);
          if(verbose) printf("k: %d - %d\n",k_min,k_max);
          if(verbose) printf("l: %d - %d\n",l_min,l_max);
          if(h_range < 0 || k_range < 0 || l_range < 0) {
              if(verbose) printf("ERROR: not enough HKL indices in %s\n",hklfilename);
              exit(9);
          }
  
          /* allocate memory for 3d arrays */
          if(verbose>6) printf("allocating %d %ld-byte double**\n",h_range+1,sizeof(double**));
          Fhkl = (double***) calloc(h_range+1,sizeof(double**));
          if(Fhkl==NULL){perror("ERROR");exit(9);};
          for (h0=0; h0<=h_range;h0++) {
              if(verbose>6) printf("allocating %d %ld-byte double*\n",k_range+1,sizeof(double*));
              Fhkl[h0] = (double**) calloc(k_range+1,sizeof(double*));
              if(Fhkl[h0]==NULL){perror("ERROR");exit(9);};
              for (k0=0; k0<=k_range;k0++) {
              if(verbose>6) printf("allocating %d %ld-byte double\n",l_range+1,sizeof(double));
                  Fhkl[h0][k0] = (double*) calloc(l_range+1,sizeof(double));
                  if(Fhkl[h0][k0]==NULL){perror("ERROR");exit(9);};
              }
          }
          if(verbose) printf("initializing to default_F = %g:\n",default_F);
          for (h0=0; h0<h_range;h0++) {
              for (k0=0; k0<k_range;k0++) {
                  for (l0=0; l0<l_range;l0++) {
                      if(verbose>9) printf("initializing %d %d %d to default_F = %g:\n",h0,k0,l0,default_F);
                      Fhkl[h0][k0][l0] = default_F;
                  }
              }
          }
          if(verbose) printf("done initializing to default_F:\n");
          if(verbose) printf("settting F000 to %g\n",F000);
          Fhkl[-h_min][-k_min][-l_min] = F000;
      }
  
      if(hklfilename != NULL)
      {
          if(verbose) printf("re-reading %s\n",hklfilename);
          while(4 == fscanf(infile,"%d%d%d%lg",&h0,&k0,&l0,&F_cell))
          {
              Fhkl[h0-h_min][k0-k_min][l0-l_min]=F_cell;
          }
          fclose(infile);
      }
  
      if(pythony_indices.size() && pythony_amplitudes.size())
      {
          if(verbose) printf("initializing Fhkl with pythony indices and amplitudes\n");
          miller_t hkl;
          for (i=0; i < pythony_indices.size(); ++i)
          {
              hkl = pythony_indices[i];
              F_cell = pythony_amplitudes[i];
              h0 = hkl[0];
              k0 = hkl[1];
              l0 = hkl[2];
              Fhkl[h0-h_min][k0-k_min][l0-l_min]=F_cell;
              if(verbose>9) printf("F %d : %d %d %d = %g\n",i,h0,k0,l0,F_cell);
          }
          if(verbose) printf("settting F000 to %g\n",F000);
          Fhkl[-h_min][-k_min][-l_min] = F000;
          if(verbose) printf("done initializing Fhkl:\n");
      }
  }
  // end of init_Fhkl
  
  
  
  
  
  /* read in or generate background profile */
  void
  nanoBragg::init_background()
  {
      int i;
      /* straighten up amorphous material properties */
      if(amorphous_sample_y<=0.0) amorphous_sample_y = beamsize;
      if(amorphous_sample_z<=0.0) amorphous_sample_z = beamsize;
      amorphous_volume = amorphous_sample_x*amorphous_sample_y*amorphous_sample_z;
      if(amorphous_density==0 && amorphous_molecules!=0 && amorphous_volume!=0 )
      {
          amorphous_density = amorphous_molecules/amorphous_volume/Avogadro*amorphous_molecular_weight;
      }
      amorphous_molecules = amorphous_volume*amorphous_density*Avogadro/amorphous_molecular_weight;
      if(verbose>1) printf("amorphous_molecules= %g in beam: %g m^3 * %g g/m^3 * %g /mol / %g g/mol\n",
               amorphous_molecules,amorphous_volume,amorphous_density,Avogadro,amorphous_molecular_weight);
  
      /* now read in amorphous material structure factors */
      stols = 0;
  
      /* option to read from a text file */
      if(stolfilename != NULL)
      {
          if(verbose) printf("reading %s\n",stolfilename);
          stols = read_text_file(stolfilename,2,&stol_of,&Fbg_of);
          if(stols == 0){
              perror("no data in input file");
              exit(9);
          }
          allocated_stols=stols;
      }
  
      /* initialize from python flex array of vec2s */
      if(pythony_stolFbg.size())
      {
          /* we will need enough space for this */
          stols = pythony_stolFbg.size();
      }
  
      if(allocated_stols != stols)
      {
          /* free any previous allocations */
          if(stol_of != NULL) {
              if(verbose>6) printf("freeing %d %ld-byte double stol_of at %p\n",stols,sizeof(double),stol_of);
              free(stol_of);
          }
          if(Fbg_of != NULL) {
              if(verbose>6) printf("freeing %d %ld-byte double Fbg_of at %p\n",stols,sizeof(double),Fbg_of);
              free(Fbg_of);
          }
          if(bin_start != NULL){
              if(verbose>6) printf("freeing %d %ld-byte double *s from bin_start at %p\n",stols,sizeof(double *),bin_start);
              free(bin_start);
          }
          if(pixels_in != NULL){
              if(verbose>6) printf("freeing %d %ld-byte unsigned ints from pixels_in at %p\n",stols,sizeof(unsigned int),pixels_in);
              free(pixels_in);
          }
  
          if(verbose>6) printf("allocating %d %ld-byte doubles for stol_of\n",stols,sizeof(double));
          stol_of = (double *) calloc(stols+10,sizeof(double));
          if(verbose>6) printf("allocating %d %ld-byte doubles for Fbg_of\n",stols,sizeof(double));
          Fbg_of  = (double *) calloc(stols+10,sizeof(double));
  
          /* allocate memory for counting how many of these get used */
          /* starting point for pixel value data for each stol-bin */
          if(verbose>6) printf("allocating %d %ld-byte double *s for bin_start\n",stols,sizeof(double *));
          bin_start = (double **) calloc(stols+2,sizeof(double *));
          /* storage for counting number of pixels in each bin */
          if(verbose>6) printf("allocating %d %ld-byte unsigned ints for pixels_in\n",stols,sizeof(unsigned int));
          pixels_in = (unsigned int *) calloc(stols,sizeof(unsigned int));
          allocated_stols = stols;
      }
  
      if(pythony_stolFbg.size())
      {
          if(verbose) printf("initializing Fbg with pythony array %d\n",stols);
          for (i=0; i < stols; ++i)
          {
              stol_of[i] = pythony_stolFbg[i][0];
              Fbg_of[i]  = pythony_stolFbg[i][1];
              if(verbose>6) printf("Fbg # %d ( stol= %g ) = %g\n",i,stol_of[i],Fbg_of[i]);
          }
          /* flag to add padding on top and bottom */
          stol_of[stols] = NAN;
          if(verbose) printf("done initializing stol_of and Fbg_of\n");
      }
  
      if(stols == 0 && amorphous_volume != 0.0)
      {
          /* do something clever here */
  
      }
  
      if(stols > 0 && isnan(stol_of[stols]))
      {
          if(verbose>9) printf("shifting F_bg vs stol array up by 2\n");
          /* clear the flag */
          stol_of[stols] = 0.0;
          /* add two values at either end for interpolation */
          stols += 4;
          Fbg_highangle = NAN;
          for(i=stols-3;i>1;--i){
              /* shift data up, and convert stol to meters while we are at it */
              stol_of[i] = stol_of[i-2] * stol_file_mult;
              Fbg_of[i]   = Fbg_of[i-2];
              if(! isnan(Fbg_of[i])) {
                  /* keep track of lowest-angle valid value */
                  Fbg_lowangle = Fbg_of[i];
                  if(isnan(Fbg_highangle)) {
                      /* also keep track of highest angle value */
                      Fbg_highangle = Fbg_of[i];
                  }
              }
              else
              {
                  /* missing values are set to default */
                  Fbg_of[i] = default_Fbg;
              }
          }
  
          /* to turn interpolation into extrapolation, bottom two values should be out at stol=negative infinity */
          stol_of[0] = -1e99;
          stol_of[1] = -1e98;
          Fbg_of[0] = Fbg_of[1] = Fbg_lowangle;
          /* to turn interpolation into extrapolation, bottom two values should be out at stol=infinity */
          stol_of[stols-2] = 1e98;
          stol_of[stols-1] = 1e99;
          Fbg_of[stols-1] = Fbg_of[stols-2] = Fbg_highangle;
      }
  }
  
  
  /* print out actual phi values in sweep */
  void
  nanoBragg::show_phisteps()
  {
      /* show phi steps with sweep over spindle axis */
      for(phi_tic = 0; phi_tic < phisteps; ++phi_tic){
          phi = phi0 + phistep*phi_tic;
          if(verbose) printf("phi%d = %g\n",phi_tic,phi*RTD);
      }
  }
  
  
  /* print out actual detector layers in sweep */
  void
  nanoBragg::show_detector_thicksteps()
  {
      /* print out detector sensor thickness with sweep over all sensor layers */
      for(thick_tic=0;thick_tic<detector_thicksteps;++thick_tic){
          if(verbose) printf("thick%d = %g um\n",thick_tic,detector_thickstep*thick_tic*1e6);
      }
  }
  
  
  /* read in or generate x-ray sources */
  void
  nanoBragg::init_sources()
  {
      /* import sources from user file */
      sources = 0;
      if(sourcefilename != NULL) {
          sources = read_text_file(sourcefilename,5,&source_X,&source_Y,&source_Z,&source_I,&source_lambda);
          if(sources == 0) {
              perror("reading source definition file");
              exit(9);
          }
          /* apply defaults to missing values */
          for(source=0;source<sources;++source){
              if(isnan(source_X[source])) {
                  source_X[source] = -source_distance*beam_vector[1];
              }
              if(isnan(source_Y[source])) {
                  source_Y[source] = -source_distance*beam_vector[2];
              }
              if(isnan(source_Z[source])) {
                  source_Z[source] = -source_distance*beam_vector[3];
              }
              if(isnan(source_I[source])) {
                  source_I[source] = 1.0;
              }
              if(isnan(source_lambda[source])) {
                  source_lambda[source] = lambda0;
              }
          }
      }
  
  
  
      if(pythony_beams.size())
      {
          if(verbose>8) printf("pythony_beams.size()= %ld\n",pythony_beams.size());
          sources = pythony_beams.size();
          if(verbose>8) printf("total sources: %d\n",sources);
          if(allocated_sources != sources && sources>0)
          {
              /* free any previous allocation */
              if(source_X != NULL) free(source_X);
              if(source_Y != NULL) free(source_Y);
              if(source_Z != NULL) free(source_Z);
              if(source_I != NULL) free(source_I);
              if(source_lambda != NULL) free(source_lambda);
  
              /* allocate enough space */
              if(verbose>6) printf("allocating space for %d sources\n",sources);
              source_X = (double *) calloc(sources+10,sizeof(double));
              source_Y = (double *) calloc(sources+10,sizeof(double));
              source_Z = (double *) calloc(sources+10,sizeof(double));
              source_I = (double *) calloc(sources+10,sizeof(double));
              source_lambda = (double *) calloc(sources+10,sizeof(double));
              allocated_sources = sources;
          }
  
          if(verbose) printf("initializing sources with pythony sources\n");
          vec3 xyz,beamdir=vec3(0,0,0),polarvec = vec3(0,0,0);
          double flux_sum = 0.0, lambda_sum = 0.0, polar_sum = 0.0, div_sum = 0.0;
          for (i=0; i < sources; ++i)
          {
              xyz = pythony_beams[i].get_sample_to_source_direction();
              source_X[i] = xyz[0];
              source_Y[i] = xyz[1];
              source_Z[i] = xyz[2];
              source_I[i] = pythony_beams[i].get_flux();
              if(isnan(source_I) || source_I[i]==0.) source_I[i] = 1.0/sources;
              flux_sum += source_I[i];
              source_lambda[i] = pythony_beams[i].get_wavelength();
              if(isnan(source_lambda[i]) || source_lambda[i]==0.) source_lambda[i] = lambda0;
              lambda_sum += source_lambda[i];
              if(verbose>8) printf("source %d :  xyz= ( %g %g %g ), I= %g lambda=%g\n",i,source_X[i],source_Y[i],source_Z[i],source_I[i],source_lambda[i]);
  
              /* average quantities that we store as one value */
              div_sum += pythony_beams[i].get_divergence();
              polar_sum += pythony_beams[i].get_polarization_fraction();
              polarvec += pythony_beams[i].get_polarization_normal();
              beamdir += xyz;
          }
          /* update averaged parameters */
          if(lambda_sum>0.0) lambda0 = lambda_sum/sources;
          polarization = polar_sum/sources;
          hdivrange=vdivrange=div_sum/sources;
          vert_vector[1] = polarvec[0];vert_vector[2] = polarvec[1];vert_vector[3] = polarvec[2];
          beam_vector[1] = beamdir[0]; beam_vector[2] = beamdir[1]; beam_vector[3] = beamdir[2];
          unitize(beam_vector,beam_vector);
          unitize(vert_vector,vert_vector);
          cross_product(beam_vector,vert_vector,polar_vector);
          unitize(polar_vector,polar_vector);
  
  
          /* take in total flux */
          if(flux_sum > 0)
          {
              flux = flux_sum;
              init_beam();
          }
          /* make sure stored source intensities are fractional */
          double norm = flux_sum/sources;
          for (i=0; i < sources && norm>0.0; ++i)
          {
              source_I[i] /= norm;
          }
  
          if(verbose) printf("done initializing sources:\n");
      }
  
  
  
      if(pythony_source_XYZ.size() || pythony_source_intensity.size() || pythony_source_lambda.size())
      {
          if(verbose>8) printf("pythony_source_XYZ.size()= %ld\n",pythony_source_XYZ.size());
          if(verbose>8) printf("pythony_source_intensity.size()= %ld\n",pythony_source_intensity.size());
          if(verbose>8) printf("pythony_source_lambda.size()= %ld\n",pythony_source_lambda.size());
          sources = pythony_source_XYZ.size();
          if(sources < pythony_source_intensity.size()) sources = pythony_source_intensity.size();
          if(sources < pythony_source_lambda.size()) sources = pythony_source_lambda.size();
          if(verbose>8) printf("total sources: %d\n",sources);
          if(allocated_sources != sources && sources>0)
          {
              /* free any previous allocation */
              if(source_X != NULL) free(source_X);
              if(source_Y != NULL) free(source_Y);
              if(source_Z != NULL) free(source_Z);
              if(source_I != NULL) free(source_I);
              if(source_lambda != NULL) free(source_lambda);
  
              /* allocate enough space */
              if(verbose>6) printf("allocating space for %d sources\n",sources);
              source_X = (double *) calloc(sources+10,sizeof(double));
              source_Y = (double *) calloc(sources+10,sizeof(double));
              source_Z = (double *) calloc(sources+10,sizeof(double));
              source_I = (double *) calloc(sources+10,sizeof(double));
              source_lambda = (double *) calloc(sources+10,sizeof(double));
              allocated_sources = sources;
          }
  
          /* make sure sizes match, or else? */
          pythony_source_XYZ.resize(sources);
          pythony_source_intensity.resize(sources);
          pythony_source_lambda.resize(sources);
  
          if(verbose) printf("initializing sources with pythony sources\n");
          vec3 xyz;
          for (i=0; i < pythony_source_XYZ.size(); ++i)
          {
              xyz = pythony_source_XYZ[i];
              source_X[i] = xyz[0];
              source_Y[i] = xyz[1];
              source_Z[i] = xyz[2];
              source_I[i] = pythony_source_intensity[i];
              if(isnan(source_I[i]) || source_I[i]==0.) source_I[i] = 1.0/sources;
              source_lambda[i] = pythony_source_lambda[i];
              if(isnan(source_lambda[i]) || source_lambda[i]==0.) source_lambda[i] = lambda0;
              if(verbose>8) printf("source %d :  xyz= ( %g %g %g ), I= %g lambda=%g\n",i,source_X[i],source_Y[i],source_Z[i],source_I[i],source_lambda[i]);
          }
          if(verbose) printf("done initializing sources:\n");
      }
  
  
      if(sources == 0)
      {
          /* generate generic list of sources */
  
          /* count divsteps sweep over solid angle of beam divergence */
          divsteps = 0;
          for(hdiv_tic=0;hdiv_tic<hdivsteps;++hdiv_tic){
              for(vdiv_tic=0;vdiv_tic<vdivsteps;++vdiv_tic){
                  hdiv = hdivstep * hdiv_tic - hdivrange/2.0 ;
                  vdiv = vdivstep * vdiv_tic - vdivrange/2.0 ;
                  /* force an elliptical divergence */
                  test = (hdiv*hdiv-hdivstep*hdivstep/4.0*(1-hdivsteps%2))/hdivrange/hdivrange ;
                  test += (vdiv*vdiv-vdivstep*vdivstep/4.0*(1-vdivsteps%2))/vdivrange/vdivrange ;
                  if( round_div && test*4.0 > 1.1) continue;
  
                  ++divsteps;
                  if(verbose) printf("divergence deviation: %g %g\n",hdiv,vdiv);
              }
          }
  
          /* print out wavelength steps with sweep over spectral dispersion */
          for(disp_tic=0;disp_tic<dispsteps;++disp_tic){
              lambda = lambda0 * ( 1.0 + dispstep * disp_tic - dispersion/2.0 ) ;
              if(verbose) printf("lambda%d = %.15g\n",disp_tic,lambda);
          }
          if(verbose) printf("dispersion= %g dispstep= %g  dispsteps= %d\n",dispersion,dispstep,dispsteps);
  
          /* free any previous allocation */
          if(source_X != NULL) free(source_X);
          if(source_Y != NULL) free(source_Y);
          if(source_Z != NULL) free(source_Z);
          if(source_I != NULL) free(source_I);
          if(source_lambda != NULL) free(source_lambda);
  
          /* allocate enough space */
          sources = divsteps*dispsteps;
          if(verbose>6) printf("allocating space for %d sources\n",sources);
          source_X = (double *) calloc(sources+10,sizeof(double));
          source_Y = (double *) calloc(sources+10,sizeof(double));
          source_Z = (double *) calloc(sources+10,sizeof(double));
          source_I = (double *) calloc(sources+10,sizeof(double));
          source_lambda = (double *) calloc(sources+10,sizeof(double));
  
          /* now actually create the source entries */
          weight = 1.0/sources;
          sources = 0;
          for(hdiv_tic=0;hdiv_tic<hdivsteps;++hdiv_tic){
              for(vdiv_tic=0;vdiv_tic<vdivsteps;++vdiv_tic){
                  hdiv = hdivstep * hdiv_tic - hdivrange/2.0 ;
                  vdiv = vdivstep * vdiv_tic - vdivrange/2.0 ;
                  /* force an elliptical divergence */
                  test = (hdiv*hdiv-hdivstep*hdivstep/4.0*(1-hdivsteps%2))/hdivrange/hdivrange ;
                  test += (vdiv*vdiv-vdivstep*vdivstep/4.0*(1-vdivsteps%2))/vdivrange/vdivrange ;
                  if( round_div && test*4.0 > 1.1) continue;
  
                  /* construct unit vector along "beam" */
                  vector[1] = -source_distance*beam_vector[1];
                  vector[2] = -source_distance*beam_vector[2];
                  vector[3] = -source_distance*beam_vector[3];
                  /* divergence is in angle space */
                  /* define "horizontal" as the E-vector of the incident beam */
                  rotate_axis(vector,newvector,polar_vector,vdiv);
                  rotate_axis(newvector,vector,vert_vector,hdiv);
  
                  /* one source at each position for each wavelength */
                  for(disp_tic=0;disp_tic<dispsteps;++disp_tic){
                      lambda = lambda0 * ( 1.0 + dispstep * disp_tic - dispersion/2.0 ) ;
  
                      source_X[sources] = vector[1];
                      source_Y[sources] = vector[2];
                      source_Z[sources] = vector[3];
                      source_I[sources] = weight;
                      source_lambda[sources] = lambda;
                      ++sources;
                  }
              }
          }
      }
      if(verbose) printf("  created a total of %d sources:\n",sources);
      for(source=0;source<sources;++source){
  
          /* retrieve stuff from cache */
          X = source_X[source];
          Y = source_Y[source];
          Z = source_Z[source];
          I = source_I[source];
          lambda = source_lambda[source];
  
          if(verbose) printf("%g %g %g   %g %.6g\n",X,Y,Z,I,lambda);
      }
  }
  // end of init_sources()
  
  
  
  
  
  /* expose internal source data */
  void
  nanoBragg::show_sources()
  {
      printf("  created a total of %d sources:\n",sources);
      for(source=0;source<sources;++source){
  
          /* retrieve stuff from cache */
          X = source_X[source];
          Y = source_Y[source];
          Z = source_Z[source];
          I = source_I[source];
          lambda = source_lambda[source];
  
          printf("%d %g %g %g   %g %.6g\n",source,X,Y,Z,I,lambda);
      }
      printf("pythony sources:\n");
      for(source=0;source<pythony_source_lambda.size();++source){
  
          /* retrieve stuff from cache */
          X = pythony_source_XYZ[source][0];
          Y = pythony_source_XYZ[source][1];
          Z = pythony_source_XYZ[source][2];
          I = pythony_source_intensity[source];
          lambda = pythony_source_lambda[source];
  
          printf("pythony: %d %g %g %g   %g %.6g\n",source,X,Y,Z,I,lambda);
      }
  }
  // end of init_sources()
  
  
  
  
  /* read in? or generate mosaic domains */
  void
  nanoBragg::init_mosaicity()
  {
      /* temporary local seed */
      long mseed;
  
      /* catch nks initialized values here? */
      if(user_mosdomains) return;
  
      /* free any previous allocations */
      if(mosaic_umats!=NULL) free(mosaic_umats);
  
      /* allocate enough space */
      if(mosaic_domains<1) mosaic_domains=1;
      if(verbose>6) printf("allocating enough space for %d mosaic domain orientation matrices\n",mosaic_domains);
      mosaic_umats = (double *) calloc(mosaic_domains+10,9*sizeof(double));
  
      /* re-initialize the RNG for the mosaic sequence */
      mseed = -labs(mosaic_seed);
  
      /* now actually create the orientation of each domain */
      for(mos_tic=0;mos_tic<mosaic_domains;++mos_tic)
      {
          if(mos_tic==0)
          {
              /* force at least one domain to be "aligned"? */
              mosaic_umats[0]=1.0;mosaic_umats[1]=0.0;mosaic_umats[2]=0.0;
              mosaic_umats[3]=0.0;mosaic_umats[4]=1.0;mosaic_umats[5]=0.0;
              mosaic_umats[6]=0.0;mosaic_umats[7]=0.0;mosaic_umats[8]=1.0;
          }
          else
          {
              mosaic_rotation_umat(mosaic_spread, mosaic_umats+9*mos_tic, &mseed);
          }
      }
  
      if(verbose) printf("  created a total of %d mosaic domains\n",mosaic_domains);
  }
  // end of init_mosaicity()
  
  
  
  
  /* print out individual mosaic domain orientations */
  void
  nanoBragg::show_mosaic_blocks()
  {
      /* assume init_mosaicity() was already run? */
      for(mos_tic=0;mos_tic<mosaic_domains;++mos_tic)
      {
  //      printf("%d diagonal %f %f %f\n",mos_tic,mosaic_umats[mos_tic*9],mosaic_umats[mos_tic*9+4],mosaic_umats[mos_tic*9+8]);
          printf("%d by: %f deg\n",mos_tic,acos((mosaic_umats[mos_tic*9]+mosaic_umats[mos_tic*9+4]+mosaic_umats[mos_tic*9+8]-1)/2)*RTD);
          umat2misset(mosaic_umats+9*mos_tic,mosaic_missets);
          printf("%d by: %f %f %f deg\n",mos_tic,mosaic_missets[1]*RTD,mosaic_missets[2]*RTD,mosaic_missets[3]*RTD);
          printf("       %f %f %f\n",*(mosaic_umats+9*mos_tic+0),*(mosaic_umats+9*mos_tic+1),*(mosaic_umats+9*mos_tic+2));
          printf("       %f %f %f\n",*(mosaic_umats+9*mos_tic+3),*(mosaic_umats+9*mos_tic+4),*(mosaic_umats+9*mos_tic+5));
          printf("       %f %f %f\n",*(mosaic_umats+9*mos_tic+6),*(mosaic_umats+9*mos_tic+7),*(mosaic_umats+9*mos_tic+8));
      }
  
      printf("  total of %d mosaic domains\n",mosaic_domains);
  }
  // end of show_mosaic_blocks()
  
  
  /* return the unitary matrices U that define the mosaic block distribution, cctbx format please */
  af::shared<mat3>
  nanoBragg::get_mosaic_blocks() {
      /* assume init_mosaicity() was already run? */
      af::shared<mat3> result;
      for(mos_tic=0;mos_tic<mosaic_domains;++mos_tic) {
          result.push_back(mat3( *(mosaic_umats+9*mos_tic+0), *(mosaic_umats+9*mos_tic+1), *(mosaic_umats+9*mos_tic+2),
                                 *(mosaic_umats+9*mos_tic+3), *(mosaic_umats+9*mos_tic+4), *(mosaic_umats+9*mos_tic+5),
                                 *(mosaic_umats+9*mos_tic+6), *(mosaic_umats+9*mos_tic+7), *(mosaic_umats+9*mos_tic+8)));
      }
      return result;
  }
  // end of get_mosaic_blocks()
  
  /* set the mosaic domains */
  void
  nanoBragg::set_mosaic_blocks(af::shared<mat3> umat_in) {
      /* free any previous allocations */
      if(mosaic_umats!=NULL) free(mosaic_umats);
  
      /* flag to avoid over-writing with automatic random domains */
      user_mosdomains=true;
  
      /* allocate enough space */
      mosaic_domains = umat_in.size();
      if(verbose>6) printf("allocating enough space for %d mosaic domain orientation matrices\n",mosaic_domains);
      mosaic_umats = (double *) calloc(mosaic_domains+10,9*sizeof(double));
  
      /* now actually import the orientation of each domain */
      for(mos_tic=0;mos_tic<mosaic_domains;++mos_tic) {
        int offset = 9 * mos_tic;
        mat3& domain = umat_in[mos_tic];
        mosaic_umats[0+offset]=domain[0];mosaic_umats[1+offset]=domain[1];mosaic_umats[2+offset]=domain[2];
        mosaic_umats[3+offset]=domain[3];mosaic_umats[4+offset]=domain[4];mosaic_umats[5+offset]=domain[5];
        mosaic_umats[6+offset]=domain[6];mosaic_umats[7+offset]=domain[7];mosaic_umats[8+offset]=domain[8];
     }
  
      if(verbose) printf("  imported a total of %d mosaic domains\n",mosaic_domains);
  }
  // end of set_mosaic_blocks()
  
  
  /* reconcile different conventions of beam center input and detector position */
  void
  nanoBragg::show_params()
  {
      printf("nanoBragg nanocrystal diffraction simulator - James Holton and Ken Frankel 8-1-17\n");
  
      printf("  %d initialized hkls (all others =%g)\n",hkls,default_F);
      printf("  ");
      if(xtal_shape == ROUND)  printf("ellipsoidal");
      if(xtal_shape == SQUARE) printf("parallelpiped");
      if(xtal_shape == GAUSS ) printf("gaussian");
      if(xtal_shape == TOPHAT) printf("tophat-spot");
      printf(" xtal: %.0fx%.0fx%.0f cells\n",Na,Nb,Nc);
      printf("Unit Cell: %g %g %g %g %g %g\n", a_A[0],b_A[0],c_A[0],alpha*RTD,beta*RTD,gamma*RTD);
      printf("Recp Cell: %g %g %g %g %g %g\n", a_star[0],b_star[0],c_star[0],alpha_star*RTD,beta_star*RTD,gamma_star*RTD);
      printf("volume = %g A^3\n",V_cell);
      printf("  amorphous part of sample is %lg m thick x %lg m high x %lg m wide, %lg g/cm^3 and %lg g/mol (%lg molecules)\n",
              amorphous_sample_x,amorphous_sample_y,amorphous_sample_z,amorphous_density/1e6,amorphous_molecular_weight,amorphous_molecules);
  
      printf("missets: %11.8f %11.8f %11.8f\n",misset[1]*RTD,misset[2]*RTD,misset[3]*RTD);
  
      /* print out the real-space matrix */
      printf("real-space cell vectors (Angstrom):\n");
      printf("     %-10s  %-10s  %-10s\n","a","b","c");
      printf("X: %11.8f %11.8f %11.8f\n",a_A[1],b_A[1],c_A[1]);
      printf("Y: %11.8f %11.8f %11.8f\n",a_A[2],b_A[2],c_A[2]);
      printf("Z: %11.8f %11.8f %11.8f\n",a_A[3],b_A[3],c_A[3]);
      printf("reciprocal-space cell vectors (Angstrom^-1):\n");
      printf("     %-10s  %-10s  %-10s\n","a_star","b_star","c_star");
      printf("X: %11.8f %11.8f %11.8f\n",a_star[1],b_star[1],c_star[1]);
      printf("Y: %11.8f %11.8f %11.8f\n",a_star[2],b_star[2],c_star[2]);
      printf("Z: %11.8f %11.8f %11.8f\n",a_star[3],b_star[3],c_star[3]);
      printf("  wave=%g meters +/- %g%% in %d steps\n",lambda0,dispersion*100,dispsteps);
      if(nopolar) { printf("  polarization effect disabled\n"); }
             else { printf("  Kahn polarization factor: %f\n",polarization); }
      if(curved_detector) printf("  curved detector: all pixels same distance from origin\n");
      if(point_pixel) printf("  pixel obliquity effect disabled\n");
      if(Fmap_pixel) printf("  Fmap mode: plotting structure factor at every pixel\n");
      printf("  incident fluence: %lg photons/m^2\n",fluence);
      printf("  distance=%lg detsize=%lgx%lg  pixel=%lg meters (%dx%d pixels)\n",distance,detsize_f,detsize_s,pixel_size,fpixels,spixels);
      printf("  Xbeam=%lg Ybeam=%lg\n",Xbeam,Ybeam);
      printf("  Fbeam=%lg Sbeam=%lg\n",Fbeam,Sbeam);
      printf("  Xclose=%lg Yclose=%lg\n",Xclose,Yclose);
      printf("  Fclose=%lg Sclose=%lg\n",Fclose,Sclose);
      printf("  DIRECTION_OF_DETECTOR_X-AXIS= %g %g %g\n",fdet_vector[1],fdet_vector[2],fdet_vector[3]);
      printf("  DIRECTION_OF_DETECTOR_Y-AXIS= %g %g %g\n",sdet_vector[1],sdet_vector[2],sdet_vector[3]);
      printf("  DIRECTION_OF_DETECTOR_Z-AXIS= %g %g %g\n",odet_vector[1],odet_vector[2],odet_vector[3]);
      printf("  INCIDENT_BEAM_DIRECTION= %g %g %g\n",beam_vector[1],beam_vector[2],beam_vector[3]);
      printf("  spindle ROTATION_AXIS= %g %g %g\n",spindle_vector[1],spindle_vector[2],spindle_vector[3]);
      cross_product(beam_vector,polar_vector,vector);
      printf("  POLARIZATION_PLANE_NORMAL= %g %g %g\n",vector[1],vector[2],vector[3]);
      printf("  dials origin= %g %g %g\n",dials_origin[1],dials_origin[2],dials_origin[3]);
      printf("  roi: %d < x < %d && %d < y < %d\n",roi_xmin,roi_xmax,roi_ymin,roi_ymax);
      printf("  hdivrange=%g hdivstep=%g  radians\n",hdivrange,hdivstep);
      printf("  vdivrange=%g vdivstep=%g  radians\n",vdivrange,vdivstep);
      printf("  %d divergence steps\n",divsteps);
      printf("  %d sources\n",sources);
      printf("  %d mosaic domains over mosaic spread of %g degrees\n",mosaic_domains,mosaic_spread*RTD);
      printf("  %d phi steps from %g to %g degrees\n",phisteps,phi0*RTD,(phi0+osc)*RTD);
      printf("  %dx%d pixel oversample steps\n",oversample,oversample);
      if(maskimage != NULL) printf("  skipping zero-flagged pixels in %s\n",maskfilename);
  //    printf("  coherent source: %d\n",coherent);
  }
  // end of show_params()
  
  
  
  
  /* add spots from nanocrystal simulation */
  void
  nanoBragg::randomize_orientation()
  {
      if(verbose>1) printf("MISSET seed %ld\n",seed);
      /* re-initialize the RNG */
      seed = -labs(seed);
      /* use spherical cap as sphere to generate random orientation in umat */
      mosaic_rotation_umat(90.0, umat, &seed);
      /* get the missetting angles, in case we want to use them again on -misset option */
      umat2misset(umat,misset);
      if(verbose) printf("random orientation misset angles: %f %f %f deg\n",misset[1]*RTD,misset[2]*RTD,misset[3]*RTD);
      /* apply this orientation shift */
      //rotate_umat(a_star,a_star,umat);
      //rotate_umat(b_star,b_star,umat);
      //rotate_umat(c_star,c_star,umat);
      /* apply below */
      misset[0] = 1.0;
      init_cell();
      return;
  }
  
  
  /* add spots from nanocrystal simulation */
  void
  nanoBragg::add_nanoBragg_spots()
  {
      max_I = 0.0;
      i = 0;
      floatimage = raw_pixels.begin();
  //    double* floatimage(raw_pixels.begin());
  //    floatimage = (double *) calloc(spixels*fpixels+10,sizeof(double));
  
      if(verbose) printf("TESTING sincg(1,1)= %f\n",sincg(1,1));
  
      /* make sure we are normalizing with the right number of sub-steps */
      steps = phisteps*mosaic_domains*oversample*oversample;
      subpixel_size = pixel_size/oversample;
  
      sum = sumsqr = 0.0;
      i = sumn = 0;
      progress_pixel = 0;
      omega_sum = 0.0;
      for(spixel=0;spixel<spixels;++spixel)
      {
          for(fpixel=0;fpixel<fpixels;++fpixel)
          {
  
              /* allow for just one part of detector to be rendered */
              if(fpixel < roi_xmin || fpixel > roi_xmax || spixel < roi_ymin || spixel > roi_ymax)
              {
                  ++i; continue;
              }
              /* allow for the use of a mask */
              if(maskimage != NULL)
              {
                  /* skip any flagged pixels in the mask */
                  if(maskimage[i] == 0)
                  {
                      ++i; continue;
                  }
              }
  
              /* reset photon count for this pixel */
              I = 0;
  
              /* loop over sub-pixels */
              for(subS=0;subS<oversample;++subS)
              {
                  for(subF=0;subF<oversample;++subF)
                  {
  
                      /* absolute mm position on detector (relative to its origin) */
                      Fdet = subpixel_size*(fpixel*oversample + subF ) + subpixel_size/2.0;
                      Sdet = subpixel_size*(spixel*oversample + subS ) + subpixel_size/2.0;
  //                  Fdet = pixel_size*fpixel;
  //                  Sdet = pixel_size*spixel;
  
                      for(thick_tic=0;thick_tic<detector_thicksteps;++thick_tic)
                      {
                          /* assume "distance" is to the front of the detector sensor layer */
                          Odet = thick_tic*detector_thickstep;
  
                          /* construct detector subpixel position in 3D space */
  //                      pixel_X = distance;
  //                      pixel_Y = Sdet-Ybeam;
  //                      pixel_Z = Fdet-Xbeam;
                          pixel_pos[1] = Fdet*fdet_vector[1]+Sdet*sdet_vector[1]+Odet*odet_vector[1]+pix0_vector[1];
                          pixel_pos[2] = Fdet*fdet_vector[2]+Sdet*sdet_vector[2]+Odet*odet_vector[2]+pix0_vector[2];
                          pixel_pos[3] = Fdet*fdet_vector[3]+Sdet*sdet_vector[3]+Odet*odet_vector[3]+pix0_vector[3];
                          pixel_pos[0] = 0.0;
                          if(curved_detector) {
                              /* construct detector pixel that is always "distance" from the sample */
                                  vector[1] = distance*beam_vector[1];
                                  vector[2]=distance*beam_vector[2] ;
                                  vector[3]=distance*beam_vector[3];
                              /* treat detector pixel coordinates as radians */
                              rotate_axis(vector,newvector,sdet_vector,pixel_pos[2]/distance);
                              rotate_axis(newvector,pixel_pos,fdet_vector,pixel_pos[3]/distance);
      //                      rotate(vector,pixel_pos,0,pixel_pos[3]/distance,pixel_pos[2]/distance);
                          }
                          /* construct the diffracted-beam unit vector to this sub-pixel */
                          airpath = unitize(pixel_pos,diffracted);
  
                          /* solid angle subtended by a pixel: (pix/airpath)^2*cos(2theta) */
                          omega_pixel = pixel_size*pixel_size/airpath/airpath*close_distance/airpath;
                          /* option to turn off obliquity effect, inverse-square-law only */
                          if(point_pixel) omega_pixel = 1.0/airpath/airpath;
                          omega_sum += omega_pixel;
  
                          /* now calculate detector thickness effects */
                          if(detector_thick > 0.0 && detector_attnlen > 0.0)
                          {
                              /* inverse of effective thickness increase */
                              parallax = dot_product(diffracted,odet_vector);
                              capture_fraction = exp(-thick_tic*detector_thickstep/detector_attnlen/parallax)
                                                -exp(-(thick_tic+1)*detector_thickstep/detector_attnlen/parallax);
                          }
                          else
                          {
                              capture_fraction = 1.0;
                          }
  
                      /* loop over sources now */
                      for(source=0;source<sources;++source){
  
                          /* retrieve stuff from cache */
                          incident[1] = -source_X[source];
                          incident[2] = -source_Y[source];
                          incident[3] = -source_Z[source];
                          lambda = source_lambda[source];
  
                          /* construct the incident beam unit vector while recovering source distance */
                          source_path = unitize(incident,incident);
  
                          /* construct the scattering vector for this pixel */
                          scattering[1] = (diffracted[1]-incident[1])/lambda;
                          scattering[2] = (diffracted[2]-incident[2])/lambda;
                          scattering[3] = (diffracted[3]-incident[3])/lambda;
  
                          /* sin(theta)/lambda is half the scattering vector length */
                          stol = 0.5*magnitude(scattering);
  
                          /* rough cut to speed things up when we aren't using whole detector */
                          if(dmin > 0.0 && stol > 0.0)
                          {
                              if(dmin > 0.5/stol)
                              {
                                  continue;
                              }
                          }
  
                          /* sweep over phi angles */
                          for(phi_tic = 0; phi_tic < phisteps; ++phi_tic)
                          {
                              phi = phi0 + phistep*phi_tic;
  
                              if( phi != 0.0 )
                              {
                                  /* rotate about spindle if neccesary */
                                  rotate_axis(a0,ap,spindle_vector,phi);
                                  rotate_axis(b0,bp,spindle_vector,phi);
                                  rotate_axis(c0,cp,spindle_vector,phi);
                              }
  
                              /* enumerate mosaic domains */
                              for(mos_tic=0;mos_tic<mosaic_domains;++mos_tic)
                              {
                                  /* apply mosaic rotation after phi rotation */
                                  if( mosaic_spread > 0.0 )
                                  {
                                      rotate_umat(ap,a,&mosaic_umats[mos_tic*9]);
                                      rotate_umat(bp,b,&mosaic_umats[mos_tic*9]);
                                      rotate_umat(cp,c,&mosaic_umats[mos_tic*9]);
                                  }
                                  else
                                  {
                                      a[1]=ap[1];a[2]=ap[2];a[3]=ap[3];
                                      b[1]=bp[1];b[2]=bp[2];b[3]=bp[3];
                                      c[1]=cp[1];c[2]=cp[2];c[3]=cp[3];
                                  }
  //                              printf("%d %f %f %f\n",mos_tic,mosaic_umats[mos_tic*9+0],mosaic_umats[mos_tic*9+1],mosaic_umats[mos_tic*9+2]);
  //                              printf("%d %f %f %f\n",mos_tic,mosaic_umats[mos_tic*9+3],mosaic_umats[mos_tic*9+4],mosaic_umats[mos_tic*9+5]);
  //                              printf("%d %f %f %f\n",mos_tic,mosaic_umats[mos_tic*9+6],mosaic_umats[mos_tic*9+7],mosaic_umats[mos_tic*9+8]);
  
                                  /* construct fractional Miller indicies */
                                  h = dot_product(a,scattering);
                                  k = dot_product(b,scattering);
                                  l = dot_product(c,scattering);
  
                                  /* round off to nearest whole index */
                                  h0 = static_cast<int>(ceil(h-0.5));
                                  k0 = static_cast<int>(ceil(k-0.5));
                                  l0 = static_cast<int>(ceil(l-0.5));
  
  
                                  /* structure factor of the lattice (paralelpiped crystal)
                                      F_latt = sin(M_PI*Na*h)*sin(M_PI*Nb*k)*sin(M_PI*Nc*l)/sin(M_PI*h)/sin(M_PI*k)/sin(M_PI*l);
                                  */
                                  F_latt = 1.0;
                                      if(xtal_shape == SQUARE)
                                  {
                                      /* xtal is a paralelpiped */
                                      if(Na>1){
                                          F_latt *= sincg(M_PI*h,Na);
                                      }
                                      if(Nb>1){
                                          F_latt *= sincg(M_PI*k,Nb);
                                      }
                                      if(Nc>1){
                                          F_latt *= sincg(M_PI*l,Nc);
                                      }
                                  }
                                      else
                                      {
                                          /* handy radius in reciprocal space, squared */
                                          hrad_sqr = (h-h0)*(h-h0)*Na*Na + (k-k0)*(k-k0)*Nb*Nb + (l-l0)*(l-l0)*Nc*Nc ;
                                      }
                                      if(xtal_shape == ROUND)
                                      {
                                          /* use sinc3 for elliptical xtal shape,
                                             correcting for sqrt of volume ratio between cube and sphere */
                                          F_latt = Na*Nb*Nc*0.723601254558268*sinc3(M_PI*sqrt( hrad_sqr * fudge ) );
                                      }
                                      if(xtal_shape == GAUSS)
                                      {
                                          /* fudge the radius so that volume and FWHM are similar to square_xtal spots */
                                          F_latt = Na*Nb*Nc*exp(-( hrad_sqr / 0.63 * fudge ));
                                      }
                                      if(xtal_shape == TOPHAT)
                                      {
                                          /* make a flat-top spot of same height and volume as square_xtal spots */
                                          F_latt = Na*Nb*Nc*(hrad_sqr*fudge < 0.3969 );
                                  }
                                  /* no need to go further if result will be zero */
                                  if(F_latt == 0.0) continue;
  
  
                                  /* find nearest point on Ewald sphere surface? */
                                  if( integral_form )
                                  {
  
                                      if( phi != 0.0 || mos_tic > 0 )
                                      {
                                          /* need to re-calculate reciprocal matrix */
  
                                          /* various cross products */
                                          cross_product(a,b,a_cross_b);
                                          cross_product(b,c,b_cross_c);
                                          cross_product(c,a,c_cross_a);
  
                                          /* new reciprocal-space cell vectors */
                                          vector_scale(b_cross_c,a_star,1e20/V_cell);
                                          vector_scale(c_cross_a,b_star,1e20/V_cell);
                                          vector_scale(a_cross_b,c_star,1e20/V_cell);
                                      }
  
                                      /* reciprocal-space coordinates of nearest relp */
                                      relp[1] = h0*a_star[1] + k0*b_star[1] + l0*c_star[1];
                                      relp[2] = h0*a_star[2] + k0*b_star[2] + l0*c_star[2];
                                      relp[3] = h0*a_star[3] + k0*b_star[3] + l0*c_star[3];
  //                                  d_star = magnitude(relp)
  
                                      /* reciprocal-space coordinates of center of Ewald sphere */
                                      Ewald0[1] = -incident[1]/lambda/1e10;
                                      Ewald0[2] = -incident[2]/lambda/1e10;
                                      Ewald0[3] = -incident[3]/lambda/1e10;
  //                                  1/lambda = magnitude(Ewald0)
  
                                       /* distance from Ewald sphere in lambda=1 units */
                                      vector[1] = relp[1]-Ewald0[1];
                                      vector[2] = relp[2]-Ewald0[2];
                                      vector[3] = relp[3]-Ewald0[3];
                                      d_r = magnitude(vector)-1.0;
  
                                      /* unit vector of diffracted ray through relp */
                                      unitize(vector,diffracted0);
  
                                      /* intersection with detector plane */
                                      xd = dot_product(fdet_vector,diffracted0);
                                      yd = dot_product(sdet_vector,diffracted0);
                                      zd = dot_product(odet_vector,diffracted0);
  
                                      /* where does the central direct-beam hit */
                                      xd0 = dot_product(fdet_vector,incident);
                                      yd0 = dot_product(sdet_vector,incident);
                                      zd0 = dot_product(odet_vector,incident);
  
                                      /* convert to mm coordinates */
                                      Fdet0 = distance*(xd/zd) + Xbeam;
                                      Sdet0 = distance*(yd/zd) + Ybeam;
  
                                      if(verbose>8) printf("integral_form: %g %g   %g %g\n",Fdet,Sdet,Fdet0,Sdet0);
                                      test = exp(-( (Fdet-Fdet0)*(Fdet-Fdet0)+(Sdet-Sdet0)*(Sdet-Sdet0) + d_r*d_r )/1e-8);
                                  } // end of integral form
  
  
                                  /* structure factor of the unit cell */
                                  if(interpolate){
                                      h0_flr = static_cast<int>(floor(h));
                                      k0_flr = static_cast<int>(floor(k));
                                      l0_flr = static_cast<int>(floor(l));
  
  
                                      if ( ((h-h_min+3)>h_range) ||
                                           (h-2<h_min)           ||
                                           ((k-k_min+3)>k_range) ||
                                           (k-2<k_min)           ||
                                           ((l-l_min+3)>l_range) ||
                                           (l-2<l_min)  ) {
                                          if(babble){
                                              babble=0;
                                              if(verbose) printf ("WARNING: out of range for three point interpolation: h,k,l,h0,k0,l0: %g,%g,%g,%d,%d,%d \n", h,k,l,h0,k0,l0);
                                              if(verbose) printf("WARNING: further warnings will not be printed! ");
                                          }
                                          F_cell = default_F;
                                          interpolate=0;
                                          continue;
                                      }
  
                                      /* integer versions of nearest HKL indicies */
                                      h_interp[0]=h0_flr-1;
                                      h_interp[1]=h0_flr;
                                      h_interp[2]=h0_flr+1;
                                      h_interp[3]=h0_flr+2;
                                      k_interp[0]=k0_flr-1;
                                      k_interp[1]=k0_flr;
                                      k_interp[2]=k0_flr+1;
                                      k_interp[3]=k0_flr+2;
                                      l_interp[0]=l0_flr-1;
                                      l_interp[1]=l0_flr;
                                      l_interp[2]=l0_flr+1;
                                      l_interp[3]=l0_flr+2;
  
                                      /* polin function needs doubles */
                                      h_interp_d[0] = (double) h_interp[0];
                                      h_interp_d[1] = (double) h_interp[1];
                                      h_interp_d[2] = (double) h_interp[2];
                                      h_interp_d[3] = (double) h_interp[3];
                                      k_interp_d[0] = (double) k_interp[0];
                                      k_interp_d[1] = (double) k_interp[1];
                                      k_interp_d[2] = (double) k_interp[2];
                                      k_interp_d[3] = (double) k_interp[3];
                                      l_interp_d[0] = (double) l_interp[0];
                                      l_interp_d[1] = (double) l_interp[1];
                                      l_interp_d[2] = (double) l_interp[2];
                                      l_interp_d[3] = (double) l_interp[3];
  
                                      /* now populate the "y" values (nearest four structure factors in each direction) */
                                      for (i1=0;i1<4;i1++) {
                                          for (i2=0;i2<4;i2++) {
                                             for (i3=0;i3<4;i3++) {
                                                    sub_Fhkl[i1][i2][i3]= Fhkl[h_interp[i1]-h_min][k_interp[i2]-k_min][l_interp[i3]-l_min];
                                             }
                                          }
                                       }
  
  
                                      /* run the tricubic polynomial interpolation */
                                      polin3(h_interp_d,k_interp_d,l_interp_d,sub_Fhkl,h,k,l,&F_cell);
                                  }
  
                                  if(! interpolate)
                                  {
                                      if ( (h0<=h_max) && (h0>=h_min) && (k0<=k_max) && (k0>=k_min) && (l0<=l_max) && (l0>=l_min)  ) {
                                          /* just take nearest-neighbor */
                                          F_cell = Fhkl[h0-h_min][k0-k_min][l0-l_min];
                                      }
                                      else
                                      {
                                          F_cell = default_F; // usually zero
                                      }
                                  }
  
                                  /* now we have the structure factor for this pixel */
  
                                  /* polarization factor */
                                  if(! nopolar){
                                      /* need to compute polarization factor */
                                      polar = polarization_factor(polarization,incident,diffracted,polar_vector);
                                  }
                                  else
                                  {
                                      polar = 1.0;
                                  }
  
                                  /* convert amplitudes into intensity (photons per steradian) */
                                      I += F_cell*F_cell*F_latt*F_latt*source_I[source]*capture_fraction*omega_pixel;
                              }
                              /* end of mosaic loop */
                          }
                          /* end of phi loop */
                      }
                      /* end of source loop */
                  }
                      /* end of detector thickness loop */
                  }
                  /* end of sub-pixel y loop */
              }
              /* end of sub-pixel x loop */
  
  
              floatimage[i] += r_e_sqr*fluence*spot_scale*polar*I/steps;
  //          floatimage[i] = test;
              if(floatimage[i] > max_I) {
                  max_I = floatimage[i];
                  max_I_x = Fdet;
                  max_I_y = Sdet;
              }
              sum += floatimage[i];
              sumsqr += floatimage[i]*floatimage[i];
              ++sumn;
  
              if( printout )
              {
                  if((fpixel==printout_fpixel && spixel==printout_spixel) || printout_fpixel < 0)
                  {
                      twotheta = atan2(sqrt(pixel_pos[2]*pixel_pos[2]+pixel_pos[3]*pixel_pos[3]),pixel_pos[1]);
                      test = sin(twotheta/2.0)/(lambda0*1e10);
                      printf("%4d %4d : stol = %g or %g\n", fpixel,spixel,stol,test);
                      printf("at %g %g %g\n", pixel_pos[1],pixel_pos[2],pixel_pos[3]);
                      printf("hkl= %f %f %f  hkl0= %d %d %d\n", h,k,l,h0,k0,l0);
                      printf(" F_cell=%g  F_latt=%g   I = %g\n", F_cell,F_latt,I);
                      printf("I/steps %15.10g\n", I/steps);
                      printf("polar   %15.10g\n", polar);
                      printf("omega   %15.10g\n", omega_pixel);
                      printf("pixel   %15.10g\n", floatimage[i]);
                      printf("real-space cell vectors (Angstrom):\n");
                      printf("     %-10s  %-10s  %-10s\n","a","b","c");
                      printf("X: %11.8f %11.8f %11.8f\n",a[1]*1e10,b[1]*1e10,c[1]*1e10);
                      printf("Y: %11.8f %11.8f %11.8f\n",a[2]*1e10,b[2]*1e10,c[2]*1e10);
                      printf("Z: %11.8f %11.8f %11.8f\n",a[3]*1e10,b[3]*1e10,c[3]*1e10);
                  }
              }
              else
              {
                  if(progress_meter && verbose && progress_pixels/100 > 0)
                  {
                      if(progress_pixel % ( progress_pixels/20 ) == 0 ||
                         ((10*progress_pixel<progress_pixels ||
                           10*progress_pixel>9*progress_pixels) &&
                          (progress_pixel % (progress_pixels/100) == 0)))
                      {
                          printf("%lu%% done\n",progress_pixel*100/progress_pixels);
                      }
                  }
                  ++progress_pixel;
              }
              ++i;
          }
      }
      if(verbose) printf("done with pixel loop\n");
  
      if(verbose) printf("solid angle subtended by detector = %g steradian ( %g%% sphere)\n",omega_sum/steps,100*omega_sum/steps/4/M_PI);
      if(verbose) printf("max_I= %g sum= %g avg= %g\n",max_I,sum,sum/sumn);
  
  }
  // end of add_nanoBragg_spots()
  
  
  
  
  /* member function to generate background from Fbg vs stol list
     arguments allow override of features that usually just slow things down,
     like oversampling pixels and multiple sources.  Providing these arguments
     does NOT change the values of the member variables */
  void
  nanoBragg::add_background( int oversample, int source )
  {
      int i;
      int source_start = 0;
      int sources = this->sources;
      max_I = 0.0;
      floatimage = raw_pixels.begin();
  //    double* floatimage(raw_pixels.begin());
  //    floatimage = (double *) calloc(spixels*fpixels+10,sizeof(double));
  
      /* allow user to override automated oversampling decision at call time with arguments */
      if(oversample<=0) oversample = this->oversample;
      if(oversample<=0) oversample = 1;
      if(source>=0) {
          /* user-specified source in the argument */
          source_start = source;
          sources = source_start +1;
      }
  
      /* make sure we are normalizing with the right number of sub-steps */
      steps = oversample*oversample;
      subpixel_size = pixel_size/oversample;
  
      /* sweep over detector */
      sum = sumsqr = 0.0;
      sumn = 0;
      progress_pixel = 0;
      omega_sum = 0.0;
      nearest = 0;
      i = 0;
      for(spixel=0;spixel<spixels;++spixel)
      {
          for(fpixel=0;fpixel<fpixels;++fpixel)
          {
              /* allow for just one part of detector to be rendered */
              if(fpixel < roi_xmin || fpixel > roi_xmax || spixel < roi_ymin || spixel > roi_ymax) {
                  invalid_pixel[i] = true;
                  ++i; continue;
              }
  
              /* reset background photon count for this pixel */
              Ibg = 0;
  
              /* loop over sub-pixels */
              for(subS=0;subS<oversample;++subS){
                  for(subF=0;subF<oversample;++subF){
  
                      /* absolute mm position on detector (relative to its origin) */
                      Fdet = subpixel_size*(fpixel*oversample + subF ) + subpixel_size/2.0;
                      Sdet = subpixel_size*(spixel*oversample + subS ) + subpixel_size/2.0;
  //                    Fdet = pixel_size*fpixel;
  //                    Sdet = pixel_size*spixel;
  
                      for(thick_tic=0;thick_tic<detector_thicksteps;++thick_tic)
                      {
                          /* assume "distance" is to the front of the detector sensor layer */
                          Odet = thick_tic*detector_thickstep;
  
                          /* construct detector pixel position in 3D space */
      //                    pixel_X = distance;
      //                    pixel_Y = Sdet-Ybeam;
      //                    pixel_Z = Fdet-Xbeam;
                          pixel_pos[1] = Fdet*fdet_vector[1]+Sdet*sdet_vector[1]+Odet*odet_vector[1]+pix0_vector[1];
                          pixel_pos[2] = Fdet*fdet_vector[2]+Sdet*sdet_vector[2]+Odet*odet_vector[2]+pix0_vector[2];
                          pixel_pos[3] = Fdet*fdet_vector[3]+Sdet*sdet_vector[3]+Odet*odet_vector[3]+pix0_vector[3];
                          pixel_pos[0] = 0.0;
                          if(curved_detector) {
                              /* construct detector pixel that is always "distance" from the sample */
                              vector[1] = distance*beam_vector[1]; vector[2]=distance*beam_vector[2] ; vector[3]=distance*beam_vector[3];
                              /* treat detector pixel coordinates as radians */
                              rotate_axis(vector,newvector,sdet_vector,pixel_pos[2]/distance);
                              rotate_axis(newvector,pixel_pos,fdet_vector,pixel_pos[3]/distance);
      //                             rotate(vector,pixel_pos,0,pixel_pos[3]/distance,pixel_pos[2]/distance);
                          }
                          /* construct the diffracted-beam unit vector to this pixel */
                          airpath = unitize(pixel_pos,diffracted);
  
                          /* solid angle subtended by a pixel: (pix/airpath)^2*cos(2theta) */
                          omega_pixel = pixel_size*pixel_size/airpath/airpath*close_distance/airpath;
                          /* option to turn off obliquity effect, inverse-square-law only */
                          if(point_pixel) omega_pixel = 1.0/airpath/airpath;
                          omega_sum += omega_pixel;
  
                          /* now calculate detector thickness effects */
                          if(detector_thick > 0.0)
                          {
                              /* inverse of effective thickness increase */
                              parallax = dot_product(diffracted,odet_vector);
                              capture_fraction = exp(-thick_tic*detector_thickstep/detector_attnlen/parallax)
                                                -exp(-(thick_tic+1)*detector_thickstep/detector_attnlen/parallax);
                          }
                          else
                          {
                              capture_fraction = 1.0;
                          }
  
                          /* loop over sources now */
                          for(source=source_start;source<sources;++source){
  
                              /* retrieve stuff from cache */
                              incident[1] = -source_X[source];
                              incident[2] = -source_Y[source];
                              incident[3] = -source_Z[source];
                              lambda = source_lambda[source];
  
                              /* construct the incident beam unit vector while recovering source distance */
                              source_path = unitize(incident,incident);
  
                              /* construct the scattering vector for this pixel */
                              scattering[1] = (diffracted[1]-incident[1])/lambda;
                              scattering[2] = (diffracted[2]-incident[2])/lambda;
                              scattering[3] = (diffracted[3]-incident[3])/lambda;
  
                              /* sin(theta)/lambda is half the scattering vector length */
                              stol = 0.5*magnitude(scattering);
  
                              /* now we need to find the nearest four "stol file" points */
                              while(stol > stol_of[nearest] && nearest <= stols){++nearest; };
                              while(stol < stol_of[nearest] && nearest >= 2){--nearest; };
  
                              /* cubic spline interpolation */
                              polint(stol_of+nearest-1, Fbg_of+nearest-1, stol, &Fbg);
  
                              /* allow negative F values to yield negative intensities */
                              sign=1.0;
                              if(Fbg<0.0) sign=-1.0;
  
                              /* now we have the structure factor for this pixel */
  
                              /* polarization factor */
                              if(! nopolar){
                                  /* need to compute polarization factor */
                                  polar = polarization_factor(polarization,incident,diffracted,polar_vector);
                              }
                              else
                              {
                                  polar = 1.0;
                              }
  
                              /* accumulate unscaled pixel intensity from this */
                              Ibg += sign*Fbg*Fbg*polar*omega_pixel*source_I[source]*capture_fraction;
                              if(verbose>7 && i==1)printf("DEBUG: Fbg= %g polar= %g omega_pixel= %g source[%d]= %g capture_fraction= %g\n",
                                                             Fbg,polar,omega_pixel,source,source_I[source],capture_fraction);
                          }
                          /* end of source loop */
                      }
                      /* end of detector thickness loop */
                  }
                  /* end of sub-pixel y loop */
              }
              /* end of sub-pixel x loop */
  
  
              /* save photons/pixel (if fluence specified), or F^2/omega if no fluence given */
              floatimage[i] += Ibg*r_e_sqr*fluence*amorphous_molecules/steps;
  
              if(verbose>7 && i==1)printf(
                "DEBUG: Ibg= %g r_e_sqr= %g fluence= %g amorphous_molecules= %g parallax= %g capfrac= %g omega= %g polar= %g steps= %d\n",
                          Ibg,r_e_sqr,fluence,amorphous_molecules,parallax,capture_fraction,omega_pixel,polar,steps);
  
              /* override: just plot interpolated structure factor at every pixel, useful for making absorption masks */
              if(Fmap_pixel) floatimage[i]= Fbg;
  
              /* keep track of basic statistics */
              if(floatimage[i] > max_I || i==0) {
                  max_I = floatimage[i];
                  max_I_x = Fdet;
                  max_I_y = Sdet;
              }
              sum += floatimage[i];
              sumsqr += floatimage[i]*floatimage[i];
              ++sumn;
  
              /* debugging infrastructure */
              if( printout )
              {
                  if((fpixel==printout_fpixel && spixel==printout_spixel) || printout_fpixel < 0)
                  {
                      twotheta = atan2(sqrt(pixel_pos[2]*pixel_pos[2]+pixel_pos[3]*pixel_pos[3]),pixel_pos[1]);
                      test = sin(twotheta/2.0)/(lambda0*1e10);
                      printf("%4d %4d : stol = %g or %g\n", fpixel,spixel,stol,test);
                      printf(" F=%g    I = %g\n", F,I);
                      printf("I/steps %15.10g\n", I/steps);
                      printf("polar   %15.10g\n", polar);
                      printf("omega   %15.10g\n", omega_pixel);
                      printf("pixel   %15.10g\n", floatimage[i]);
                  }
              }
              else
              {
                  if(progress_meter && progress_pixels/100 > 0)
                  {
                      if(progress_pixel % ( progress_pixels/20 ) == 0 ||
                         ((10*progress_pixel<progress_pixels ||
                           10*progress_pixel>9*progress_pixels) &&
                          (progress_pixel % (progress_pixels/100) == 0)))
                      {
                          printf("%lu%% done\n",progress_pixel*100/progress_pixels);
                      }
                  }
                  ++progress_pixel;
              }
              /* end progress meter stuff */
  
              /* never ever forget to increment this */
              ++i;
          }
          /* end fpixel loop */
      }
      /* end spixel loop */
  
      if(verbose) printf("\nsolid angle subtended by detector = %g steradian ( %g%% sphere)\n",omega_sum/steps,100*omega_sum/steps/4/M_PI);
      if(verbose) printf("max_I= %g @ ( %g, %g) sum= %g avg= %g\n",max_I,max_I_x,max_I_y,sum,sum/sumn);
  
  }
  // end of add_background()
  
  
  
  
  
  
  /* function for extracting stol-vs-Fbg data from an existing image */
  void
  nanoBragg::extract_background(int source)
  {
      /* make sure we override k from hkl */
      int i,k;
  
      nearest = 0;
      max_I = 0.0;
      floatimage = raw_pixels.begin();
  //    double* floatimage(raw_pixels.begin());
  //    floatimage = (double *) calloc(spixels*fpixels+10,sizeof(double));
  
      /* override internal default if argument was given */
      if(source<0) {
          /* user did not specify source in the argument, use the first one */
          source = 0;
      }
      /* no oversampling, this is an extraction */
      steps = 1;
  
      /* sweep over detector */
      sum = sumsqr = 0.0;
      i = 0;
      progress_pixel = 0;
      valid_pixels = 0;
      omega_sum = 0.0;
      for(spixel=0;spixel<spixels;++spixel)
      {
          for(fpixel=0;fpixel<fpixels;++fpixel)
          {
  
              /* allow for just one part of detector to be rendered */
              if(fpixel < roi_xmin || fpixel > roi_xmax || spixel < roi_ymin || spixel > roi_ymax) {
                  invalid_pixel[i] = true;
                  ++i; continue;
              }
  
              /* is the pixel valid on the input image? */
              /* skip over any invalid values */
              for(k=1;k<=ignore_values;++k)
              {
                  if(floatimage[i]==ignore_value[k]){
                      invalid_pixel[i] = true;
                  }
              }
  //            if(invalid_pixel[i])
  //            {
  //                ++i;
  //                continue;
  //            }
  
              /* still need polar, omega and stol, but cannot deconvolute sub-pixels and source size,
                 so assume oversample=1, neutronium detector, and take the first source */
              thick_tic = 0;
              //source = 0;
  
              /* absolute mm position on detector (relative to its origin) */
              //Fdet = subpixel_size*(fpixel*oversample + subF ) + subpixel_size/2.0;
              //Sdet = subpixel_size*(spixel*oversample + subS ) + subpixel_size/2.0;
              Fdet = pixel_size*fpixel;
              Sdet = pixel_size*spixel;
  
              /* assume "distance" is to the front of the detector sensor layer */
              //Odet = thick_tic*detector_thickstep;
              Odet = 0.;
  
              /* construct detector pixel position in 3D space */
              //                    pixel_X = distance;
              //                    pixel_Y = Sdet-Ybeam;
              //                    pixel_Z = Fdet-Xbeam;
              pixel_pos[1] = Fdet*fdet_vector[1]+Sdet*sdet_vector[1]+Odet*odet_vector[1]+pix0_vector[1];
              pixel_pos[2] = Fdet*fdet_vector[2]+Sdet*sdet_vector[2]+Odet*odet_vector[2]+pix0_vector[2];
              pixel_pos[3] = Fdet*fdet_vector[3]+Sdet*sdet_vector[3]+Odet*odet_vector[3]+pix0_vector[3];
              pixel_pos[0] = 0.0;
              if(curved_detector) {
                  /* construct detector pixel that is always "distance" from the sample */
                  vector[1] = distance*beam_vector[1];
                  vector[2] = distance*beam_vector[2] ;
                  vector[3] = distance*beam_vector[3];
                  /* treat detector pixel coordinates as radians */
                  rotate_axis(vector,newvector,sdet_vector,pixel_pos[2]/distance);
                  rotate_axis(newvector,pixel_pos,fdet_vector,pixel_pos[3]/distance);
                  //rotate(vector,pixel_pos,0,pixel_pos[3]/distance,pixel_pos[2]/distance);
              }
              /* construct the diffracted-beam unit vector to this pixel */
              airpath = unitize(pixel_pos,diffracted);
  
              /* solid angle subtended by a pixel: (pix/airpath)^2*cos(2theta) */
              omega_pixel = pixel_size*pixel_size/airpath/airpath*close_distance/airpath;
              /* option to turn off obliquity effect, inverse-square-law only */
              if(point_pixel) omega_pixel = 1.0/airpath/airpath;
              omega_sum += omega_pixel;
  
              /* retrieve stuff from cache */
              incident[1] = -source_X[source];
              incident[2] = -source_Y[source];
              incident[3] = -source_Z[source];
              lambda = source_lambda[source];
  
              /* construct the incident beam unit vector while recovering source distance */
              source_path = unitize(incident,incident);
  
              /* construct the scattering vector for this pixel */
              scattering[1] = (diffracted[1]-incident[1])/lambda;
              scattering[2] = (diffracted[2]-incident[2])/lambda;
              scattering[3] = (diffracted[3]-incident[3])/lambda;
  
              /* sin(theta)/lambda is half the scattering vector length */
              stol = 0.5*magnitude(scattering);
  
              /* now we need to find the nearest four "stol file" points */
              while(stol > stol_of[nearest] && nearest <= stols){++nearest; };
              while(stol < stol_of[nearest] && nearest >= 2){--nearest; };
  
              /* polarization factor */
              if(! nopolar){
                  /* need to compute polarization factor */
                  polar = polarization_factor(polarization,incident,diffracted,polar_vector);
              }
              else
              {
                  polar = 1.0;
              }
  
              /* now we have everything we need to transform pixel intensity back to a structure factor */
              /* note that for real image data we want to subtract ADC offset, but not for half-done simulated data */
              deviate=floatimage[i]-adc_offset;
              /* for negative intensities, make the structure factor negative too */
              sign = 1.0;
              if(deviate<0.0) sign = -1.0;
              deviate = fabs(deviate);
              pixel_F = sign*sqrt(deviate/polar/omega_pixel/fluence/r_e_sqr/amorphous_molecules);
              /* maintain F and stol images */
              stolimage[i] = stol/stol_file_mult;
              Fimage[i] = pixel_F;
              bin = 0;
              if(! invalid_pixel[i])
              {
                  /* figure out which stol bin this pixel belongs to.  invalid pixels are in bin=0 */
                  bin = nearest;
                  /* move to next bin if it is more than halfway there */
                  if(stol > (stol_of[bin]+stol_of[bin+1])/2.0) ++bin;
                  /* only pixels with bin assignments are valid */
                  ++valid_pixels;
              }
              ++pixels_in[bin];
              bin_of[i]=bin;
  
              /* now maybe do some stats? */
              if(floatimage[i] > max_I) {
                  max_I = floatimage[i];
                  max_I_x = Fdet;
                  max_I_y = Sdet;
              }
              sum += floatimage[i];
              sumsqr += floatimage[i]*floatimage[i];
              ++n;
  
              /* debugging option: print out particular pixel, or progress meter */
              if( printout )
              {
                  if((fpixel==printout_fpixel && spixel==printout_spixel) || printout_fpixel < 0)
                  {
                      twotheta = atan2(sqrt(pixel_pos[2]*pixel_pos[2]+pixel_pos[3]*pixel_pos[3]),pixel_pos[1]);
                      test = sin(twotheta/2.0)/(lambda0*1e10);
                          printf("%4d %4d : stol = %g or %g\n", fpixel,spixel,stol,test);
                          printf(" F=%g    I = %g\n", F,I);
                          printf("I/steps %15.10g\n", I/steps);
                          printf("polar   %15.10g\n", polar);
                          printf("omega   %15.10g\n", omega_pixel);
                          printf("pixel   %15.10g\n", floatimage[i]);
                  }
              }
              else
              {
                  if(progress_meter && progress_pixels/100 > 0)
                  {
                      if(progress_pixel % ( progress_pixels/20 ) == 0 ||
                         ((10*progress_pixel<progress_pixels ||
                           10*progress_pixel>9*progress_pixels) &&
                          (progress_pixel % (progress_pixels/100) == 0)))
                      {
                          printf("%lu%% done\n",progress_pixel*100/progress_pixels);
                      }
                  }
                  ++progress_pixel;
              }
  
              ++i;
          }
          /* end fpixel loop */
      }
      /* end spixel loop */
  
      /* now we need to organize Fpixel data into bins */
  
      /* set up pointers with enough space after each of them (2*n for median/mad filter) */
      bin_start[0]= (double *) calloc(2*pixels+10*stols,sizeof(double));
      ++bin_start[0];
      for(bin=1;bin<stols-1;++bin)
      {
          /* each array must have 2*n values in it */
          /* we counted number of pixels_in each bin in above loop */
          bin_start[bin]=bin_start[bin-1]+2*pixels_in[bin-1]+2;
      }
  
      /* populate each bin with appropriate pixel values */
      for(i=0;i<pixels;++i)
      {
          /* recover which bin this pixel is in */
          bin = bin_of[i];
          /* copy F value into start of bin data */
          *bin_start[bin] = Fimage[i];
          /* increment the pointer to the next value for this bin */
          ++bin_start[bin];
          /* yes, I know.  The bin_start now points to the end, we will reset the starting points in the next loop */
      }
  
  
      /* now we go through bins, rejecting outliers and taking mean of whats left */
      for(bin=2;bin<stols-2;++bin)
      {
          /* correct pointer drift in last loop, bin_start will now point to the start of its data */
          bin_start[bin] -= pixels_in[bin];
  
          /* recover cached value of stol for this bin */
          stol = stol_of[bin];
          /* this function looks at "input_n" elements, starting at 1 */
          median   = fmedian_with_rejection(pixels_in[bin],bin_start[bin]-1,6.0,&mad,&n);
          /* now that outliers are rejected, take the mean of what is left */
          avg_arej = fmean_with_rejection(n,bin_start[bin],6.0,&rmsd_arej,&n);
          if(n>100)
          {
  //          fprintf(outfile,"%g %g %g  %d\n",stol/stol_file_mult,median,mad,n);
  //          fprintf(outfile,"%g %g %g  %d\n",stol/stol_file_mult,avg_arej,rmsd_arej,n);
              if(verbose) printf("bin= %d stol= %g avg= %g rmsd= %g  n= %d\n",bin,stol/stol_file_mult,avg_arej,rmsd_arej,n);
              /* now populate the Fbg array.  This should automagically get passed back to Python? */
              Fbg_of[bin]=avg_arej;
          }
          else
          {
              if(verbose) printf("WARNING: not enough pixels in bin= %d n=%d stol= %g median= %g avg_arej= %g\n",bin,n,stol/stol_file_mult,median,avg_arej);
          }
      }
  
      if(verbose) printf("done with radial median filter\n");
  
      /* now we need to copy back into pythony array */
  }
  // end of extract_background()
  
  
  
  
  
  
  /* function for applying the PSF, copies over raw pixels with blurred version of itself */
  void
  nanoBragg::apply_psf(shapetype psf_type, double fwhm_pixels, int user_psf_radius)
  {
      max_I=0.0;
      double* inimage(raw_pixels.begin());
      double *outimage=NULL;
      double *kernel;
      int x0,y0,x,y,dx,dy;
      double g,rsq;
      double photon_noise,lost_photons=0.0,total_lost_photons=0.0;
      int maxwidth,kernel_size,psf_radius;
      int i,j,k;
      double photonloss_factor = 10.0; // inverse of maximum tolerable number of lost photons
  
  
      /* take the member value for PSF radius if it is set and nothing was in the function call */
      if(user_psf_radius <= 0 && this->psf_radius > 0) user_psf_radius = this->psf_radius;
      /* find a fwhm for the PSF in pixel units */
      if(fwhm_pixels <= 0.0) fwhm_pixels = this->psf_fwhm/this->pixel_size;
      /* update the members, for posterity */
      this->psf_fwhm = fwhm_pixels * this->pixel_size;
      this->psf_radius = user_psf_radius;
  
      if(verbose>7) printf("apply_psf(): user_psf_radius = %d\n",user_psf_radius);
      if(verbose>7) printf("apply_psf(): updated psf_fwhm = %g  pixel_size= %g\n",psf_fwhm,pixel_size);
  
      /* convert fwhm to "g" distance : fwhm = sqrt((2**(2./3)-1))/2*g */
      g = fwhm_pixels * 0.652383013252053;
  
      if(psf_type != GAUSS && psf_type != FIBER)
      {
          if(verbose) printf("ERROR: unknown PSF type\n");
          return;
      }
  
      pixels = fpixels*spixels;
      if(pixels == 0)
      {
          if(verbose) printf("ERROR: apply_psf image has zero size\n");
          return;
      }
  
      if(fwhm_pixels <= 0.0)
      {
          if(verbose) printf("WARNING: apply_psf function has zero size\n");
          return;
      }
  
      /* start with a clean slate */
      if(outimage!=NULL) free(outimage);
      outimage = (double *) calloc(pixels+10,sizeof(double));
  
      psf_radius = user_psf_radius;
      if(psf_radius <= 0)
      {
          /* auto-select radius */
  
          /* preliminary stats */
          max_I = 0.0;
          for(i=0;i<pixels;++i)
          {
              /* optionally scale the input file */
              if(max_I < inimage[i]) max_I = inimage[i];
          }
          if(verbose) printf("  maximum input photon/pixel: %g\n",max_I);
  
          if(max_I<=0.0)
          {
              /* nothing to blur */
              if(verbose) printf("WARNING: no photons, PSF skipped\n");
              return;
          }
  
          /* at what level will an error in intensity be lost? */
          photon_noise = sqrt(max_I);
          lost_photons = photon_noise/photonloss_factor;
          if(verbose) printf("apply_psf() predicting %g lost photons\n",lost_photons);
  
          if(psf_type == GAUSS)
          {
              /* calculate the radius beyond which only 0.5 photons will fall */
              psf_radius = 1+ceil( sqrt(-log(lost_photons/max_I)/log(4.0)/2.0)*fwhm_pixels );
              if(verbose) printf("  auto-selected psf_radius = %d x %d pixels for rendering kernel\n",psf_radius,psf_radius);
          }
          if(psf_type == FIBER)
          {
              /* calculate the radius r beyond which only 0.5 photons will fall */
              /* r = sqrt((g*(max_I/0.5))**2-g**2)
                   ~ 2*g*max_I */
              psf_radius = 1+ceil( g*(max_I/lost_photons)  );
              if(verbose) printf("  auto-selected psf_radius = %d x %d pixels for rendering kernel\n",psf_radius,psf_radius);
          }
          if(psf_radius == 0) psf_radius = 1;
      }
      /* limit psf kernel to be no bigger than 4x the input image */
      maxwidth = fpixels;
      if(spixels > maxwidth) maxwidth = spixels;
      if(psf_radius > maxwidth) psf_radius = maxwidth;
      kernel_size = 2*psf_radius+1;
      if(verbose>6) printf("apply_psf() kernel_size= %d\n",kernel_size);
  
      /* now alocate enough space to store the PSF kernel image */
      kernel = (double *) calloc(kernel_size*kernel_size,sizeof(double));
      if(kernel == NULL)
      {
          perror("apply_psf: could not allocate memory for PSF kernel");
          exit(9);
      }
  
      /* cache the PSF in an array */
      for(dy=-psf_radius;dy<=psf_radius;++dy)
      {
          for(dx=-psf_radius;dx<=psf_radius;++dx)
          {
              rsq = dx*dx+dy*dy;
              if(rsq > psf_radius*psf_radius) continue;
  
              /* this could be more efficient */
              k = kernel_size*(kernel_size/2+dy)+kernel_size/2+dx;
  
  
              if( psf_type == GAUSS ) {
                  kernel[k] = integrate_gauss_over_pixel(dx,dy,fwhm_pixels,1.0);
              }
              if( psf_type == FIBER ) {
                  kernel[k] = integrate_fiber_over_pixel(dx,dy,g,1.0);
              }
          }
      }
  
      /* implement PSF  */
      double sum_in = 0.0, sum_out = 0.0;
      sumn = 0;
      for(i=0;i<pixels;++i)
      {
          x0 = i%fpixels;
          y0 = (i-x0)/fpixels;
  
          /* skip if there is nothing to add */
          if(inimage[i] <= 0.0) continue;
  
          sum_in += inimage[i]; ++sumn;
  
          if(user_psf_radius != 0)
          {
              psf_radius = user_psf_radius;
          }
          else
          {
              /* at what level will an error in intensity be lost? */
              photon_noise = sqrt(inimage[i]);
              lost_photons = photon_noise/photonloss_factor;
  
              if(psf_type == GAUSS)
              {
                  /* calculate the radius beyond which only 0.5 photons will fall
                     r = sqrt(-log(lost_photons/total_photons)/log(4)/2)*fwhm */
                  psf_radius = 1+ceil( sqrt(-log(lost_photons/inimage[i])/log(16.0))*fwhm_pixels );
  //              printf("  auto-selected psf_radius = %d pixels\n",psf_radius);
              }
              if(psf_type == FIBER)
              {
                  /* calculate the radius beyond which only 0.5 photons will fall
                     r = sqrt((g*(total_photons/lost_photons))**2-g**2)
                       ~ g*total_photons/lost_photons */
                  psf_radius = 1+ceil( g*(inimage[i]/lost_photons)  );
  //              printf("  (%d,%d) auto-selected psf_radius = %d pixels\n",x0,y0,psf_radius);
              }
          }
          if(psf_radius == 0) psf_radius = 1;
          /* limit psf kernel to be no bigger than 4x the input image */
          maxwidth = fpixels;
          if(spixels > maxwidth) maxwidth = spixels;
          if(psf_radius > maxwidth) psf_radius = maxwidth;
  
          /* given the radius, how many photons will escape? */
          if(psf_type == GAUSS)
          {
              /* r = sqrt(-log(lost_photons/total_photons)/log(16))*fwhm */
              /* lost_photons = total_photons*exp(-log(16)*(r^2/fwhm^2)) */
              rsq = psf_radius;
              rsq = rsq/fwhm_pixels;
              rsq = rsq*rsq;
              lost_photons = inimage[i]*exp(-log(16.0)*rsq);
          }
          if(psf_type == FIBER)
          {
              /* r ~ g*total_photons/lost_photons
                 normalized integral from r=inf to "r" :  g/sqrt(g**2+r**2) */
              lost_photons = inimage[i]*g/sqrt(g*g+psf_radius*psf_radius);
          }
          /* accumulate this so we can add it to the whole image */
          total_lost_photons += lost_photons;
  
          for(dx=-psf_radius;dx<=psf_radius;++dx)
          {
              for(dy=-psf_radius;dy<=psf_radius;++dy)
              {
                  /* this could be more efficient */
                  k = kernel_size*(kernel_size/2+dy)+kernel_size/2+dx;
                  if(kernel[k] == 0.0) continue;
  
                  rsq = dx*dx+dy*dy;
                  if(rsq > psf_radius*psf_radius) continue;
                  x = x0+dx;
                  y = y0+dy;
                  if(x<0 || x>fpixels) continue;
                  if(y<0 || y>spixels) continue;
  
                  /* index into output array */
                  j = y*fpixels+x;
                  /* do not wander off the output array */
                  if(j<0 || j > pixels) continue;
  
                  outimage[j] += inimage[i]*kernel[k];
              }
          }
      }
      /* now we have some lost photons, add them back "everywhere" */
      lost_photons = total_lost_photons/pixels;
      if(verbose) printf("adding back %g lost photons\n",total_lost_photons);
      for(i=0;i<pixels;++i)
      {
          sum_out += outimage[i];
          outimage[i] += lost_photons;
          if(verbose>7 && i==pixels/2) printf("apply_psf() pixel=%d in= %g out= %g \n",i,inimage[i],outimage[i]);
      }
      if(verbose>7) printf("apply_psf() sum_in=  %g\n",sum_in);
      if(verbose>7) printf("apply_psf() sum_out= %g\n",sum_out);
      if(verbose>7) printf("apply_psf() sum_out= %g (after correction)\n",sum_out+total_lost_photons);
  
      /* don't need kernel anymore. */
      free(kernel);
  
      i=pixels/2;
      if(verbose>7) printf("apply_psf() pixel=%d in= %g out= %g \n",i,inimage[i],outimage[i]);
  
      /* and now.  No idea how to exchange buffers without confusing Python, so lets just copy it back */
      memcpy(inimage,outimage,pixels*sizeof(double));
  
      free(outimage);
      return;
  }
  // end of apply_psf()
  
  
  
  
  // function to add different types of noise to the image
  void
  nanoBragg::add_noise()
  {
      int i = 0;
      long cseed;
  
      double expected_photons,observed_photons,adu;
      /* refer to raw_pixels pixel data */
      floatimage = raw_pixels.begin();
  
      /* don't bother with this loop if calibration is perfect
       NOTE: applying calibration before Poisson noise simulates loss of photons before the detector
       NOTE: applying calibration after Poisson noise simulates systematics in read-out electronics
       here we do the latter */
  
      /* re-start the RNG */
      seed = -labs(seed);
  
      if(verbose) printf("applying calibration at %g%%, flicker noise at %g%%\n",calibration_noise*100.,flicker_noise*100.);
      sum = max_I = 0.0;
      i = sumn = 0;
      for(spixel=0;spixel<spixels;++spixel)
      {
          for(fpixel=0;fpixel<fpixels;++fpixel)
          {
  
              /* allow for just one part of detector to be rendered */
              if(fpixel < roi_xmin || fpixel > roi_xmax || spixel < roi_ymin || spixel > roi_ymax)
              {
                  ++i; continue;
              }
              /* allow for the use of a mask */
              if(maskimage != NULL)
              {
                  /* skip any flagged pixels in the mask */
                  if(maskimage[i] == 0)
                  {
                      ++i; continue;
                  }
              }
  
              /* take input image to be ideal photons/pixel */
              expected_photons = floatimage[i];
  
              /* negative photons should be taken as invalid? */
              if(expected_photons < 0.0)
              {
                  ++i; continue;
              }
  
              /* simulate 1/f noise in source */
              if(flicker_noise > 0.0){
                  expected_photons *= ( 1.0 + flicker_noise * gaussdev( &seed ) );
              }
              /* simulate photon-counting error */
              observed_photons = poidev( expected_photons, &seed );
  
              /* now we overwrite the flex array, it is now observed, rather than expected photons */
              floatimage[i] = observed_photons;
  
              /* accumulate number of photons, and keep track of max */
              if(floatimage[i] > max_I) {
                  max_I = floatimage[i];
                  max_I_x = fpixel;
                  max_I_y = spixel;
              }
              sum += observed_photons;
              ++sumn;
  
              ++i;
          }
      }
      if(verbose) printf("%.0f photons generated on noise image, max= %f at ( %.0f, %.0f )\n",sum,max_I,max_I_x,max_I_y);
  
  
  
      if(calibration_noise > 0.0)
      {
          /* calibration is same from shot to shot, so use well-known seed */
          cseed = -labs(calib_seed);
          sum = max_I = 0.0;
          i = sumn = 0;
          for(spixel=0;spixel<spixels;++spixel)
          {
              for(fpixel=0;fpixel<fpixels;++fpixel)
              {
                  /* allow for just one part of detector to be rendered */
                  if(fpixel < roi_xmin || fpixel > roi_xmax || spixel < roi_ymin || spixel > roi_ymax)
                  {
                      ++i; continue;
                  }
                  /* allow for the use of a mask */
                  if(maskimage != NULL)
                  {
                      /* skip any flagged pixels in the mask */
                      if(maskimage[i] == 0)
                      {
                          ++i; continue;
                      }
                  }
  
                  /* calibration is same from shot to shot, but varies from pixel to pixel */
                  floatimage[i] *= ( 1.0 + calibration_noise * gaussdev( &cseed ) );
  
                  /* accumulate number of photons, and keep track of max */
                  if(floatimage[i] > max_I) {
                      max_I = floatimage[i];
                      max_I_x = fpixel;
                      max_I_y = spixel;
                  }
                  sum += floatimage[i];
                  ++sumn;
  
                  ++i;
              }
          }
      }
      if(verbose) printf("%.0f photons after calibration error, max= %f at ( %.0f, %.0f )\n",sum,max_I,max_I_x,max_I_y);
  
  
      /* now would be a good time to implement PSF?  before we add read-out noise */
  
      /* now that we have photon count at each point, implement any PSF */
      if(psf_type != UNKNOWN && psf_fwhm > 0.0)
      {
          /* report on sum before the PSF is applied */
          if(verbose) printf("%.0f photons on noise image before PSF\n",sum);
          /* start with a clean slate */
          if(verbose) printf("  applying PSF width = %g um\n",psf_fwhm*1e6);
  
          apply_psf(psf_type, psf_fwhm/pixel_size, 0);
  
          /* the flex array is now the blurred version of itself, ready for read-out noise */
      }
  
  
      if(verbose) printf("adu = quantum_gain= %g * observed_photons + offset= %g + readout_noise= %g\n",quantum_gain,adc_offset,readout_noise);
      sum = max_I = 0.0;
      i = sumn = 0;
      for(spixel=0;spixel<spixels;++spixel)
      {
          for(fpixel=0;fpixel<fpixels;++fpixel)
          {
              /* allow for just one part of detector to be rendered */
              if(fpixel < roi_xmin || fpixel > roi_xmax || spixel < roi_ymin || spixel > roi_ymax)
              {
                  ++i; continue;
              }
              /* allow for the use of a mask */
              if(maskimage != NULL)
              {
                  /* skip any flagged pixels in the mask */
                  if(maskimage[i] == 0)
                  {
                      ++i; continue;
                  }
              }
  
                  /* convert photon signal to pixel units */
                  adu = floatimage[i]*quantum_gain + adc_offset;
  
                  /* readout noise is in pixel units (adu) */
                  if(readout_noise > 0.0){
                      adu += readout_noise * gaussdev( &seed );
              }
  
              /* once again, overwriting flex array, this time in ADU units */
              floatimage[i] = adu;
  
              if(adu > max_I) {
                  max_I = adu;
                  max_I_x = fpixel;
                  max_I_y = spixel;
              }
              sum += adu;
              ++sumn;
              ++i;
          }
      }
      if(verbose) printf("%.0f net adu generated on final image, max= %f at ( %.0f, %.0f )\n",sum-adc_offset*sumn,max_I,max_I_x,max_I_y);
  }
  // end of add_noise()
  
  void
  nanoBragg::to_smv_format_streambuf(boost_adaptbx::python::streambuf & output,
      double intfile_scale, int const&debug_x, int const& debug_y) const {
      boost_adaptbx::python::streambuf::ostream os(output);
      const double* floatimage = raw_pixels.begin();
      double max_value = (double)std::numeric_limits<unsigned short int>::max();
      double saturation = floor(max_value - 1 );
      /* output as ints */
  
      unsigned short int intimage;
      double max_I = this-> max_I;
      double max_I_x = this-> max_I_x;
      double max_I_y = this-> max_I_y;
      if(intfile_scale <= 0.0){
          /* need to auto-scale */
          int i=0;
          for(int spixel=0;spixel<spixels;++spixel)
          {
              for(int fpixel=0;fpixel<fpixels;++fpixel)
              {
                  if(fpixel==debug_x && spixel==debug_y) printf("DEBUG: pixel # %d at (%d,%d) has value %g\n",i,fpixel,spixel,floatimage[i]);
                  if(i==0 || max_I < floatimage[i])
                  {
                      max_I = floatimage[i];
                      max_I_x = fpixel;
                      max_I_y = spixel;
                  }
                  ++i;
              }
          }
          if(verbose) printf("providing default scaling: max_I = %g @ (%g %g)\n",max_I,max_I_x,max_I_y);
          intfile_scale = 1.0;
          if(max_I>0.0) intfile_scale = 55000.0/(max_I);
      }
      if(verbose) printf("scaling data by: intfile_scale = %g\n",intfile_scale);
  
      double sum = 0.0;
      max_I = 0.0;
      int i = 0;
      for(int spixel=0;spixel<spixels;++spixel)
      {
          for(int fpixel=0;fpixel<fpixels;++fpixel)
          {
              /* no noise, just use intfile_scale */
  
              intimage = (unsigned short int) (std::min(saturation, floatimage[i]*intfile_scale ));
              os.write((char *) &intimage, sizeof(unsigned short int));
  
              if(verbose>90) printf("DEBUG #%d %g -> %g -> %d\n",i,floatimage[i],floatimage[i]*intfile_scale,intimage);
  
              if((double) intimage > max_I || i==0) {
                  max_I = (double) intimage;
                  max_I_x = fpixel;
                  max_I_y = spixel;
              }
              if(fpixel==debug_x && spixel==debug_y) printf("DEBUG: pixel # %d at (%d,%d) has int value %d\n",i,fpixel,spixel,intimage);
  
              sum += intimage;
              ++i;
          }
      }
      // os << "Hello world";
  
  }
  
  void
  nanoBragg::to_smv_format(
      std::string const& fileout, double intfile_scale, int debug_x, int debug_y)
  {
      pixels = spixels * fpixels;
      floatimage = raw_pixels.begin();
      FILE* outfile;
      double max_value = (double)std::numeric_limits<unsigned short int>::max();
      double saturation = floor(max_value - 1 );
      static const char *byte_order = "little_endian";
  
      /* output as ints */
      af::versa<unsigned short int, af::c_grid<2> > intimage_v(
        af::c_grid<2> (spixels,fpixels));
      unsigned short int * intimage = intimage_v.begin();
  
      if(intfile_scale <= 0.0)
      {
          /* need to auto-scale */
          i=0;
          for(spixel=0;spixel<spixels;++spixel)
          {
              for(fpixel=0;fpixel<fpixels;++fpixel)
              {
                  if(fpixel==debug_x && spixel==debug_y) printf("DEBUG: pixel # %d at (%d,%d) has value %g\n",i,fpixel,spixel,floatimage[i]);
                  if(i==0 || max_I < floatimage[i])
                  {
                      max_I = floatimage[i];
                      max_I_x = fpixel;
                      max_I_y = spixel;
                  }
                  ++i;
              }
          }
          if(verbose) printf("providing default scaling: max_I = %g @ (%g %g)\n",max_I,max_I_x,max_I_y);
          intfile_scale = 1.0;
          if(max_I>0.0) intfile_scale = 55000.0/(max_I);
      }
      if(verbose) printf("scaling data by: intfile_scale = %g\n",intfile_scale);
  
      sum = max_I = 0.0;
      int i = 0;
      for(spixel=0;spixel<spixels;++spixel)
      {
          for(fpixel=0;fpixel<fpixels;++fpixel)
          {
  
              /* no noise, just use intfile_scale */
              intimage[i] = (unsigned short int) (std::min(saturation,
                             floatimage[i]*intfile_scale ));
  
              if(verbose>90) printf("DEBUG #%d %g -> %g -> %d\n",i,floatimage[i],floatimage[i]*intfile_scale,intimage[i]);
  
              if((double) intimage[i] > max_I || i==0) {
                  max_I = (double) intimage[i];
                  max_I_x = fpixel;
                  max_I_y = spixel;
              }
              if(fpixel==debug_x && spixel==debug_y) printf("DEBUG: pixel # %d at (%d,%d) has int value %d\n",i,fpixel,spixel,intimage[i]);
  
              sum += intimage[i];
              ++i;
          }
      }
      if (verbose){
        printf("writing %s as %d-byte integers; sum= %g max= %g @ ( %g %g)\n",fileout.c_str(),
              (int)sizeof(unsigned short int),sum,max_I,max_I_x,max_I_y);
      }
      outfile = fopen(fileout.c_str(),"wb");
      fprintf(outfile,"{\nHEADER_BYTES=1024;\nDIM=2;\nBYTE_ORDER=%s;\nTYPE=unsigned_short;\n",byte_order);
      fprintf(outfile,"SIZE1=%d;\nSIZE2=%d;\nPIXEL_SIZE=%g;\nDISTANCE=%g;\n",fpixels,spixels,pixel_size*1000.0,distance*1000.0);
      fprintf(outfile,"WAVELENGTH=%g;\n",lambda0*1e10);
      fprintf(outfile,"BEAM_CENTER_X=%g;\nBEAM_CENTER_Y=%g;\n",Xbeam*1000.0,Ybeam*1000);
      fprintf(outfile,"ADXV_CENTER_X=%g;\nADXV_CENTER_Y=%g;\n",Fbeam*1000.0,(detsize_s-Sbeam)*1000);
      fprintf(outfile,"MOSFLM_CENTER_X=%g;\nMOSFLM_CENTER_Y=%g;\n",(Sbeam-0.5*pixel_size)*1000.0,(Fbeam-0.5*pixel_size)*1000);
      fprintf(outfile,"DENZO_X_BEAM=%g;\nDENZO_Y_BEAM=%g;\n",(Sbeam+0.0*pixel_size)*1000.0,(Fbeam+0.0*pixel_size)*1000);
      fprintf(outfile,"DIALS_ORIGIN=%g,%g,%g\n",dials_origin[1],dials_origin[2],dials_origin[3]);
      fprintf(outfile,"XDS_ORGX=%g;\nXDS_ORGY=%g;\n",ORGX,ORGY);
      fprintf(outfile,"CLOSE_DISTANCE=%g;\n",close_distance*1000.0);
      fprintf(outfile,"PHI=%g;\nOSC_START=%g;\nOSC_RANGE=%g;\n",phi0*RTD,phi0*RTD,osc*RTD);
      fprintf(outfile,"TIME=%g;\n",exposure);
      fprintf(outfile,"TWOTHETA=%g;\n",detector_twotheta*RTD);
      fprintf(outfile,"DETECTOR_SN=000;\n");
      fprintf(outfile,"ADC_OFFSET=%g;\n",adc_offset);
      fprintf(outfile,"BEAMLINE=fake;\n");
      fprintf(outfile,"}\f");
      while ( ftell(outfile) < 1024 ){ fprintf(outfile," "); };
      fwrite(intimage,sizeof(unsigned short int),pixels,outfile);
  
      fclose(outfile);
  
      return;
  }
  
  
  
  
  /* Fourier transform of a grating */
  double sincg(double x,double N) {
      if(x==0.0) return N;
  
      return sin(x*N)/sin(x);
  }
  
  /* Fourier transform of a sphere */
  double sinc3(double x) {
      if(x==0.0) return 1.0;
  
      return 3.0*(sin(x)/x-cos(x))/(x*x);
  }
  
  double sinc_conv_sinc3(double x) {
      if(x==0.0) return 1.0;
  
      return 3.0*(sin(x)-x*cos(x))/(x*x*x);
  }
  
  
  double *rotate(double *v, double *newv, double phix, double phiy, double phiz) {
  
      double rxx,rxy,rxz,ryx,ryy,ryz,rzx,rzy,rzz;
      double new_x,new_y,new_z,rotated_x,rotated_y,rotated_z;
  
      new_x=v[1];
      new_y=v[2];
      new_z=v[3];
  
      if(phix != 0){
          /* rotate around x axis */
          //rxx= 1;         rxy= 0;         rxz= 0;
          ryx= 0;         ryy= cos(phix); ryz=-sin(phix);
          rzx= 0;         rzy= sin(phix); rzz= cos(phix);
  
          rotated_x = new_x;
          rotated_y = new_y*ryy + new_z*ryz;
          rotated_z = new_y*rzy + new_z*rzz;
          new_x = rotated_x; new_y = rotated_y; new_z = rotated_z;
      }
  
      if(phiy != 0) {
          /* rotate around y axis */
          rxx= cos(phiy); rxy= 0;         rxz= sin(phiy);
          //ryx= 0;         ryy= 1;         ryz= 0;
          rzx=-sin(phiy); rzy= 0;         rzz= cos(phiy);
  
          rotated_x = new_x*rxx + new_y*rxy + new_z*rxz;
          rotated_y = new_y;
          rotated_z = new_x*rzx + new_y*rzy + new_z*rzz;
          new_x = rotated_x; new_y = rotated_y; new_z = rotated_z;
      }
  
      if(phiz != 0){
          /* rotate around z axis */
          rxx= cos(phiz); rxy=-sin(phiz); rxz= 0;
          ryx= sin(phiz); ryy= cos(phiz); ryz= 0;
          //rzx= 0;         rzy= 0;         rzz= 1;
  
          rotated_x = new_x*rxx + new_y*rxy ;
          rotated_y = new_x*ryx + new_y*ryy;
          rotated_z = new_z;
          new_x = rotated_x; new_y = rotated_y; new_z = rotated_z;
      }
  
      newv[1]=new_x;
      newv[2]=new_y;
      newv[3]=new_z;
  
      return newv;
  }
  
  
  
  /* rotate a point about a unit vector axis */
  double *rotate_axis(double *v, double *newv, double *axis, double phi) {
  
      double sinphi = sin(phi);
      double cosphi = cos(phi);
      double dot = (axis[1]*v[1]+axis[2]*v[2]+axis[3]*v[3])*(1.0-cosphi);
      double temp[4];
  
      temp[1] = axis[1]*dot+v[1]*cosphi+(-axis[3]*v[2]+axis[2]*v[3])*sinphi;
      temp[2] = axis[2]*dot+v[2]*cosphi+(+axis[3]*v[1]-axis[1]*v[3])*sinphi;
      temp[3] = axis[3]*dot+v[3]*cosphi+(-axis[2]*v[1]+axis[1]*v[2])*sinphi;
      newv[1]=temp[1]; newv[2]=temp[2]; newv[3]=temp[3];
  
      return newv;
  }
  
  
  
  /* rotate a vector using a 9-element unitary matrix */
  double *rotate_umat(double *v, double *newv, double umat[9]) {
  
      double uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz;
  
      /* for convenience, assign matrix x-y coordinate */
      uxx = umat[0];
      uxy = umat[1];
      uxz = umat[2];
      uyx = umat[3];
      uyy = umat[4];
      uyz = umat[5];
      uzx = umat[6];
      uzy = umat[7];
      uzz = umat[8];
  
      /* rotate the vector (x=1,y=2,z=3) */
      newv[1] = uxx*v[1] + uxy*v[2] + uxz*v[3];
      newv[2] = uyx*v[1] + uyy*v[2] + uyz*v[3];
      newv[3] = uzx*v[1] + uzy*v[2] + uzz*v[3];
  
      return newv;
  }
  
  
  
  
  /* returns a 9-element unitary matrix for a random isotropic rotation on a spherical cap of diameter "mosaicity" */
  /* mosaic = 90 deg is a full sphere */
  double *mosaic_rotation_umat(double mosaicity, double umat[9], long *seed)
  {
  //    double ran1(long *idum);
      double r1,r2,r3,xyrad,rot;
      double v1,v2,v3;
      double t1,t2,t3,t6,t7,t8,t9,t11,t12,t15,t19,t20,t24;
      double uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz;
  
      /* make three random uniform deviates on [-1:1] */
      r1= (double) 2.0*ran1(seed)-1.0;
      r2= (double) 2.0*ran1(seed)-1.0;
      r3= (double) 2.0*ran1(seed)-1.0;
  
      xyrad = sqrt(1.0-r2*r2);
      rot = mosaicity*powf((1.0-r3*r3),(1.0/3.0));
  
      v1 = xyrad*sin(M_PI*r1);
      v2 = xyrad*cos(M_PI*r1);
      v3 = r2;
  
      /* commence incomprehensible quaternion calculation */
      t1 =  cos(rot);
      t2 =  1.0 - t1;
      t3 =  v1*v1;
      t6 =  t2*v1;
      t7 =  t6*v2;
      t8 =  sin(rot);
      t9 =  t8*v3;
      t11 = t6*v3;
      t12 = t8*v2;
      t15 = v2*v2;
      t19 = t2*v2*v3;
      t20 = t8*v1;
      t24 = v3*v3;
  
      /* populate the unitary rotation matrix */
      umat[0] = uxx = t1 + t2*t3;
      umat[1] = uxy = t7 - t9;
      umat[2] = uxz = t11 + t12;
      umat[3] = uyx = t7 + t9;
      umat[4] = uyy = t1 + t2*t15;
      umat[5] = uyz = t19 - t20;
      umat[6] = uzx = t11 - t12;
      umat[7] = uzy = t19 + t20;
      umat[8] = uzz = t1 + t2*t24;
  
      /* return pointer to the provided array, in case that is useful */
      return umat;
  }
  
  /* convert a unitary rotation matrix into misseting angles
     rotx roty rotz are returned as missets[1] missets[2] missets[3] */
  double *umat2misset(double umat[9],double *missets)
  {
      double uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz;
      double m,mx,my,mz;
      double xcy_x,xcy_y,xcy_z;
      double ycz_x,ycz_y,ycz_z;
      double zcx_x,zcx_y,zcx_z;
      double rotx,roty,rotz;
  
      uxx=umat[0];uxy=umat[1];uxz=umat[2];
      uyx=umat[3];uyy=umat[4];uyz=umat[5];
      uzx=umat[6];uzy=umat[7];uzz=umat[8];
  
      /* or transpose? */
  //    uxx=umat[1];uyx=umat[2];uzx=umat[3];
  //    uxy=umat[4];uyy=umat[5];uzy=umat[6];
  //    uxz=umat[7];uyz=umat[8];uzz=umat[9];
  
      /* make sure it is unitary */
      mx = sqrt(uxx*uxx+uxy*uxy+uxz*uxz);
      my = sqrt(uyx*uyx+uyy*uyy+uyz*uyz);
      mz = sqrt(uzx*uzx+uzy*uzy+uzz*uzz);
      if(mx>0){uxx/=mx;uxy/=mx;uxz/=mx;};
      if(my>0){uyx/=my;uyy/=my;uyz/=my;};
      if(mz>0){uzx/=mz;uzy/=mz;uzz/=mz;};
  
      if(mx>=0 && my<=0 && mz<=0)
      {
          uyx=0;uyy=1;uyz=0;
          uzx=0;uzy=0;uzz=1;
      }
      if(mx<=0 && my>=0 && mz<=0)
      {
          uxx=1;uxy=0;uxz=0;
          uzx=0;uzy=0;uzz=1;
      }
      if(mx<=0 && my<=0 && mz>=0)
      {
          uxx=1;uxy=0;uxz=0;
          uyx=0;uyy=1;uyz=0;
      }
  
      /* cross products to check normality */
      xcy_x = uxy*uyz - uxz*uyy;
      xcy_y = uxz*uyx - uxx*uyz;
      xcy_z = uxx*uyy - uxy*uyx;
      m=sqrt(xcy_x*xcy_x+xcy_y*xcy_y+xcy_z*xcy_z);
      if(m>0){xcy_x/=m;xcy_y/=m;xcy_z/=m;};
  
      ycz_x = uyy*uzz - uyz*uzy;
      ycz_y = uyz*uzx - uyx*uzz;
      ycz_z = uyx*uzy - uyy*uzx;
      m=sqrt(ycz_x*ycz_x+ycz_y*ycz_y+ycz_z*ycz_z);
      if(m>0){ycz_x/=m;ycz_y/=m;ycz_z/=m;};
  
      zcx_x = uzy*uxz - uzz*uxy;
      zcx_y = uzz*uxx - uzx*uxz;
      zcx_z = uzx*uxy - uzy*uxx;
      m=sqrt(zcx_x*zcx_x+zcx_y*zcx_y+zcx_z*zcx_z);
      if(m>0){zcx_x/=m;zcx_y/=m;zcx_z/=m;};
  
      /* substitute any empty vectors for cross-product of other two */
      if(mx<=0){uxx=ycz_x;uxy=ycz_y;uxz=ycz_z;};
      if(my<=0){uyx=zcx_x;uyy=zcx_y;uyz=zcx_z;};
      if(mz<=0){uzx=xcy_x;uzy=xcy_y;uzz=xcy_z;};
  
  
      /* cross products to check normality */
      xcy_x = uxy*uyz - uxz*uyy;
      xcy_y = uxz*uyx - uxx*uyz;
      xcy_z = uxx*uyy - uxy*uyx;
      m=sqrt(xcy_x*xcy_x+xcy_y*xcy_y+xcy_z*xcy_z);
      if(m>0){xcy_x/=m;xcy_y/=m;xcy_z/=m;}
  
      ycz_x = uyy*uzz - uyz*uzy;
      ycz_y = uyz*uzx - uyx*uzz;
      ycz_z = uyx*uzy - uyy*uzx;
      m=sqrt(ycz_x*ycz_x+ycz_y*ycz_y+ycz_z*ycz_z);
      if(m>0){ycz_x/=m;ycz_y/=m;ycz_z/=m;};
  
      zcx_x = uzy*uxz - uzz*uxy;
      zcx_y = uzz*uxx - uzx*uxz;
      zcx_z = uzx*uxy - uzy*uxx;
      m=sqrt(zcx_x*zcx_x+zcx_y*zcx_y+zcx_z*zcx_z);
      if(m>0){zcx_x/=m;zcx_y/=m;zcx_z/=m;};
  
  
      /* substitute any empty vectors for cross-product of other two */
      if(mx<=0){uxx=ycz_x;uxy=ycz_y;uxz=ycz_z;};
      if(my<=0){uyx=zcx_x;uyy=zcx_y;uyz=zcx_z;};
      if(mz<=0){uzx=xcy_x;uzy=xcy_y;uzz=xcy_z;};
  
  
  
      /* make sure it is unitary */
      mx = sqrt(uxx*uxx+uxy*uxy+uxz*uxz);
      my = sqrt(uyx*uyx+uyy*uyy+uyz*uyz);
      mz = sqrt(uzx*uzx+uzy*uzy+uzz*uzz);
      if(mx>0){uxx/=mx;uxy/=mx;uxz/=mx;};
      if(my>0){uyx/=my;uyy/=my;uyz/=my;};
      if(mz>0){uzx/=mz;uzy/=mz;uzz/=mz;};
  
      /* see if its really orthonormal? */
  
      if(uzx*uzx < 1.0)
      {
          rotx = atan2(uzy,uzz);
          roty = atan2(-uzx,sqrt(uzy*uzy+uzz*uzz));
          rotz = atan2(uyx,uxx);
      }
      else
      {
          rotx = atan2(1.0,1.0)*4.0;
          roty = atan2(1.0,1.0)*2.0;
          rotz = atan2(uxy,-uyy);
      }
  
      missets[1] = rotx;
      missets[2] = roty;
      missets[3] = rotz;
      return missets;
  }
  
  
  
  /* random number generators */
  
  /* Poisson deviate given expectation value of photon count (xm) */
  double poidev(double xm, long *idum)
  {
  //    double gammln(double xx);
  //    double ran1(long *idum);
      /* oldm is a flag for whether xm has changed since last call */
      static double sq,alxm,g,oldm=(-1.0);
      double em,t,y;
  
      /* routine below locks up for > 1e6 photons? */
      if (xm > 1.0e6) {
          return xm+sqrt(xm)*gaussdev(idum);
      }
  
      if (xm < 12.0) {
          /* use direct method: simulate exponential delays between events */
          if(xm != oldm) {
              /* xm is new, compute the exponential */
              oldm=xm;
              g=exp(-xm);
          }
          /* adding exponential deviates is equivalent to multiplying uniform deviates */
          /* final comparison is to the pre-computed exponential */
          em = -1;
          t = 1.0;
          do {
              ++em;
              t *= ran1(idum);
          } while (t > g);
      } else {
          /* Use rejection method */
          if(xm != oldm) {
              /* xm has changed, pre-compute a few things... */
              oldm=xm;
              sq=sqrt(2.0*xm);
              alxm=log(xm);
              g=xm*alxm-gammln(xm+1.0);
          }
          do {
              do {
                  /* y is a deviate from a lorentzian comparison function */
                  y=tan(M_PI*ran1(idum));
                  /* shift and scale */
                  em=sq*y+xm;
              } while (em < 0.0);         /* there are no negative Poisson deviates */
              /* round off to nearest integer */
              em=floor(em);
              /* ratio of Poisson distribution to comparison function */
              /* scale it back by 0.9 to make sure t is never > 1.0 */
              t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
          } while (ran1(idum) > t);
      }
  
      return em;
  }
  
  
  /* return gaussian deviate with rms=1 and FWHM = 2/sqrt(log(2)) */
  double gaussdev(long *idum)
  {
  //    double ran1(long *idum);
      static int iset=0;
      static double gset;
      double fac,rsq,v1,v2;
  
      if (iset == 0) {
          /* no extra deviats handy ... */
  
          /* so pick two uniform deviates on [-1:1] */
          do {
              v1=2.0*ran1(idum)-1.0;
              v2=2.0*ran1(idum)-1.0;
              rsq=v1*v1+v2*v2;
          } while (rsq >= 1.0 || rsq == 0);
          /* restrained to the unit circle */
  
          /* apply Box-Muller transformation to convert to a normal deviate */
          fac=sqrt(-2.0*log(rsq)/rsq);
          gset=v1*fac;
          iset=1;         /* we now have a spare deviate */
          return v2*fac;
      } else {
          /* there is an extra deviate in gset */
          iset=0;
          return gset;
      }
  }
  
  
  /* generate Lorentzian deviate with FWHM = 2 */
  double lorentzdev(long *seed) {
  //    double ran1(long *idum);
  
      return tan(M_PI*(ran1(seed)-0.5));
  }
  
  /* return triangular deviate with FWHM = 1 */
  double triangledev(long *seed) {
      double ran1(long *idum);
      double value;
  
      value = ran1(seed);
      if(value > 0.5){
          value = sqrt(2*(value-0.5))-1;
      }else{
          value = 1-sqrt(2*value);
      }
  
      return value;
  }
  
  
  
  double expdev(long *idum)
  {
      double dum;
  
      do
      dum=ran1(idum);
      while( dum == 0.0);
      return -log(dum);
  }
  
  
  
  /* ln of the gamma function */
  double gammln(double xx)
  {
      double x,y,tmp,ser;
      static double cof[6]={76.18009172947146,-86.50532032941677,
      24.01409824083091,-1.231739572450155,
      0.1208650973866179e-2,-0.5395239384953e-5};
      int j;
  
      y=x=xx;
      tmp=x+5.5;
      tmp -= (x+0.5)*log(tmp);
      ser = 1.000000000190015;
      for(j=0;j<=5;++j) ser += cof[j]/++y;
  
      return -tmp+log(2.5066282746310005*ser/x);
  }
  
  
  
  
  
  /* returns a uniform random deviate between 0 and 1 */
  #define IA 16807
  #define IM 2147483647
  #define AM (1.0/IM)
  #define IQ 127773
  #define IR 2836
  #define NTAB 32
  #define NDIV (1+(IM-1)/NTAB)
  #define EPS 1.2e-7
  #define RNMX (1.0-EPS)
  
  double ran1(long *idum)
  {
      int j;
      long k;
      static long iy=0;
      static long iv[NTAB];
      double temp;
  
      if (*idum <= 0 || !iy) {
          /* first time around.  don't want idum=0 */
          if(-(*idum) < 1) *idum=1;
          else *idum = -(*idum);
  
          /* load the shuffle table */
          for(j=NTAB+7;j>=0;j--) {
              k=(*idum)/IQ;
              *idum=IA*(*idum-k*IQ)-IR*k;
              if(*idum < 0) *idum += IM;
              if(j < NTAB) iv[j] = *idum;
          }
          iy=iv[0];
      }
      /* always start here after initializing */
      k=(*idum)/IQ;
      *idum=IA*(*idum-k*IQ)-IR*k;
      if (*idum < 0) *idum += IM;
      j=iy/NDIV;
      iy=iv[j];
      iv[j] = *idum;
      if((temp=AM*iy) > RNMX) return RNMX;
      else return temp;
  }
  
  
  void polint(double *xa, double *ya, double x, double *y)
  {
          double x0,x1,x2,x3;
          x0 = (x-xa[1])*(x-xa[2])*(x-xa[3])*ya[0]/((xa[0]-xa[1])*(xa[0]-xa[2])*(xa[0]-xa[3]));
          x1 = (x-xa[0])*(x-xa[2])*(x-xa[3])*ya[1]/((xa[1]-xa[0])*(xa[1]-xa[2])*(xa[1]-xa[3]));
          x2 = (x-xa[0])*(x-xa[1])*(x-xa[3])*ya[2]/((xa[2]-xa[0])*(xa[2]-xa[1])*(xa[2]-xa[3]));
          x3 = (x-xa[0])*(x-xa[1])*(x-xa[2])*ya[3]/((xa[3]-xa[0])*(xa[3]-xa[1])*(xa[3]-xa[2]));
          *y = x0+x1+x2+x3;
  }
  
  
  
  void polin2(double *x1a, double *x2a, double **ya, double x1, double x2, double *y)
  {
          void polint(double *xa, double *ya, double x, double *y);
          int j;
          double ymtmp[4];
          for (j=1;j<=4;j++) {
                  polint(x2a,ya[j-1],x2,&ymtmp[j-1]);
          }
          polint(x1a,ymtmp,x1,y);
  }
  
  
  void polin3(double *x1a, double *x2a, double *x3a, double ***ya, double x1,
          double x2, double x3, double *y)
  {
          void polint(double *xa, double ya[], double x, double *y);
          void polin2(double *x1a, double *x2a, double **ya, double x1,double x2, double *y);
          void polin1(double *x1a, double *ya, double x1, double *y);
          int j;
          double ymtmp[4];
  
          for (j=1;j<=4;j++) {
              polin2(x2a,x3a,&ya[j-1][0],x2,x3,&ymtmp[j-1]);
          }
          polint(x1a,ymtmp,x1,y);
  }
  
  
  /* FWHM = integral = 1 */
  double ngauss2D(double x,double y)
  {
      return log(16.0)/M_PI*exp(-log(16.0)*(x*x+y*y));
  }
  double ngauss2Dinteg(double x,double y)
  {
      return 0.125*(erf(2.0*x*sqrt(log(2.0)))*erf(y*sqrt(log(16.0)))*sqrt(log(16.0)/log(2.0)));
  }
  
  
  
  
  
  
  
  
  
  /* read in multi-column text file to list of double arrays */
  /* provide address of undeclared arrays on command line */
  size_t read_text_file(char *filename, size_t nargs, ... )
  {
      /* maximum of 10240-character lines? */
      char text[10240];
      char *token;
      const char delimiters[] = " \t,;:!";
      const char numberstuf[] = "0123456789-+.EGeg";
  
      unsigned long line,lines;
      unsigned long i,j;
      double value;
      double *data;
      double **pointer;
      va_list arglist;
      FILE *infile = NULL;
  
      infile = fopen(filename,"r");
      if(infile == NULL) {
          perror("fopen()");
          return 0;
      }
      lines=0;
      while ( fgets ( text, sizeof text, infile ) != NULL ) {
          token = text;
          token += strspn(token,delimiters);
          if(strcmp(token,"\n")==0) {
              //printf("blank\n");
              continue;
          }
          ++lines;
      }
      rewind(infile);
  
      /* allocate memory for arrays */
      va_start( arglist, nargs);
      for(i=0;i<nargs;++i){
          /* allocate the array */
          data = (double*) malloc((lines+10)*sizeof(double));
          /* initialize with missing number flags */
          for(j=0;j<lines+10;++j) {
              data[j] = NAN;
          }
          /* get argument (pointer to pointer) */
          pointer = va_arg(arglist, double **);
          /* change the value of what the arg points to */
          *pointer = data;
          /* now the pointer provided as an argument points to
          something */
      }
      va_end(arglist);
  
      line = 0;
      while ( fgets ( text, sizeof text, infile ) != NULL ) { /* read a line */
  
          token = text;
          token += strspn(token,delimiters);
          if(strcmp(token,"\n")==0) {
              //printf("blank\n");
              continue;
          }
          i=0;
          va_start( arglist, nargs);
          do
          {
              value=atof(token);
              /* get argument */
              pointer = va_arg(arglist, double **);
              /* retrieve data array's address */
              data = *pointer;
              data[line] = value;
  
              token += strspn(token,numberstuf);
              if (strcmp(token,"\n")==0) continue;
              token += strcspn(token,delimiters);
              token += strspn(token,delimiters);
              if (strcmp(token,"\n")==0) continue;
  
              ++i;
              if(i>=nargs) {
                  break;
              }
          }
          while (strcmp(token,"\n")!=0) ;
          va_end(arglist);
  
  //      printf("initializing:");
  //        va_start( arglist, nargs);
  //        for(i=0;i<nargs;++i){
  //          pointer = va_arg(arglist, double **);
  //          data = *pointer;
  //          printf(" %g",data[line]);
  //        }
  //        va_end(arglist);
  //      printf("\n");
  
          ++line;
      }
      fclose(infile);
  
      return lines;
  }
  
  
  /* measure magnitude of provided vector */
  double magnitude(double *vector) {
  
      /* measure the magnitude */
      vector[0] = sqrt(vector[1]*vector[1]+vector[2]*vector[2]+vector[3]*vector[3]);
  
      return vector[0];
  }
  
  /* make provided vector a unit vector */
  double unitize(double *vector, double *new_unit_vector) {
      double mag;
  
      /* measure the magnitude */
      mag = magnitude(vector);
  
      if(mag != 0.0){
          /* normalize it */
          new_unit_vector[1]=vector[1]/mag;
          new_unit_vector[2]=vector[2]/mag;
          new_unit_vector[3]=vector[3]/mag;
      }
      else
      {
          /* can't normalize, report zero vector */
          new_unit_vector[0] = 0.0;
          new_unit_vector[1] = 0.0;
          new_unit_vector[2] = 0.0;
          new_unit_vector[3] = 0.0;
      }
      return mag;
  }
  
  /* scale magnitude of provided vector */
  double vector_scale(double *vector, double *new_vector, double scale) {
  
      new_vector[1] = scale*vector[1];
      new_vector[2] = scale*vector[2];
      new_vector[3] = scale*vector[3];
  
      return magnitude(new_vector);
  }
  
  /* enforce magnitude of provided vector */
  double vector_rescale(double *vector, double *new_vector, double new_magnitude) {
      double oldmag;
  
      oldmag = magnitude(vector);
      if(oldmag <= 0.0) oldmag = 1.0;
      new_vector[1] = new_magnitude/oldmag*vector[1];
      new_vector[2] = new_magnitude/oldmag*vector[2];
      new_vector[3] = new_magnitude/oldmag*vector[3];
  
      return magnitude(new_vector);
  }
  
  /* difference between two given vectors */
  double vector_diff(double *vector, double *origin_vector, double *new_vector) {
  
      new_vector[1] = vector[1]-origin_vector[1];
      new_vector[2] = vector[2]-origin_vector[2];
      new_vector[3] = vector[3]-origin_vector[3];
      return magnitude(new_vector);
  }
  
  
  /* vector cross product where vector magnitude is 0th element */
  double *cross_product(double *x, double *y, double *z) {
      z[1] = x[2]*y[3] - x[3]*y[2];
      z[2] = x[3]*y[1] - x[1]*y[3];
      z[3] = x[1]*y[2] - x[2]*y[1];
      z[0] = 0.0;
  
      return z;
  }
  /* vector inner product where vector magnitude is 0th element */
  double dot_product(double *x, double *y) {
      return x[1]*y[1]+x[2]*y[2]+x[3]*y[3];
  }
  
  
  /* polarization factor */
  double polarization_factor(double kahn_factor, double *incident, double *diffracted, double *axis)
  {
      double cos2theta,cos2theta_sqr,sin2theta_sqr;
      double psi=0.0;
      double E_in[4];
      double B_in[4];
      double E_out[4];
      double B_out[4];
  
      unitize(incident,incident);
      unitize(diffracted,diffracted);
      unitize(axis,axis);
  
      /* component of diffracted unit vector along incident beam unit vector */
      cos2theta = dot_product(incident,diffracted);
      cos2theta_sqr = cos2theta*cos2theta;
      sin2theta_sqr = 1-cos2theta_sqr;
  
      if(kahn_factor != 0.0){
          /* tricky bit here is deciding which direciton the E-vector lies in for each source
             here we assume it is closest to the "axis" defined above */
  
          /* cross product to get "vertical" axis that is orthogonal to the cannonical "polarization" */
          cross_product(axis,incident,B_in);
          /* make it a unit vector */
          unitize(B_in,B_in);
  
          /* cross product with incident beam to get E-vector direction */
          cross_product(incident,B_in,E_in);
          /* make it a unit vector */
          unitize(E_in,E_in);
  
          /* get components of diffracted ray projected onto the E-B plane */
          E_out[0] = dot_product(diffracted,E_in);
          B_out[0] = dot_product(diffracted,B_in);
  
          /* compute the angle of the diffracted ray projected onto the incident E-B plane */
          psi = -atan2(B_out[0],E_out[0]);
      }
  
      /* correction for polarized incident beam */
      return 0.5*(1.0 + cos2theta_sqr - kahn_factor*cos(2*psi)*sin2theta_sqr);
  }
  
  
  
  
  
  /* 2D Gaussian integral=1 */
  double ngauss2D(double x, double y, double fwhm)
  {
      return log(16.0)/M_PI*fwhm*fwhm*exp(-log(16.0)*((x*x+y*y)/(fwhm*fwhm) ));
  }
  
  /* integral of Gaussian fwhm=1 integral=1 */
  double ngauss2D_integ(double x, double y)
  {
      return 0.125*(erf(2.0*x*sqrt(log(2.0)))*erf(y*sqrt(log(16.0)))*sqrt(log(16.0)/log(2.0)));
  }
  
  /* unit volume integrated over a pixel, fwhm = 1 */
  double ngauss2D_pixel(double x,double y,double pix)
  {
      return ngauss2D_integ(x+pix/2.,y+pix/2.)-ngauss2D_integ(x+pix/2.,y-pix/2.)-ngauss2D_integ(x-pix/2.,y+pix/2.)+ngauss2D_integ(x-pix/2.,y-pix/2.);
  }
  
  double integrate_gauss_over_pixel(double x, double y, double fwhm, double pix)
  {
      return ngauss2D_pixel(x/fwhm,y/fwhm,pix/fwhm);
  }
  
  
  double fiber2D(double x,double y,double g)
  {
      /* g/(2*pi)*(g**2+x**2+y**2)**(-3/2) */
      double temp;
      temp = sqrt(g*g+x*x+y*y);
      if(temp <= 0.0) return 0.0;
      return g/2.0/M_PI/temp/temp/temp;
  }
  double fiber2D_integ(double x,double y,double g)
  {
      return atan((x*y)/(g*sqrt(g*g + x*x + y*y)))/2.0/M_PI;
  }
  double fiber2D_pixel(double x,double y,double g,double pix)
  {
       return fiber2D_integ(x+pix/2.,y+pix/2.,g)-fiber2D_integ(x+pix/2.,y-pix/2.,g)-fiber2D_integ(x-pix/2.,y+pix/2.,g)+fiber2D_integ(x-pix/2.,y-pix/2.,g);
  }
  double integrate_fiber_over_pixel(double x, double y, double g, double pix)
  {
      return fiber2D_pixel(x,y,g,pix);
  }
  
  
  
  #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
  double fmedian(unsigned int n, double arr[])
  {
      unsigned int i,j,k,l,ir,mid;
      double a,temp;
  
      l=1;
      ir=n;
      k=(n+1)/2;
  //printf("n=%d; k=%d\n",n,k);
  
  //for(i=1;i<=n;++i) printf("arr[%d]=%f\n",i,arr[i]);
  
      for(;;)
      {
          if(ir <= l+1)
          {
              if(ir == l+1 && arr[ir] < arr[l])
              {
                  SWAP(arr[l],arr[ir]);
              }
  //for(i=1;i<=n;++i) printf("arr[%d]=%f\n",i,arr[i]);
              return arr[k];
          } else {
              mid=(l+ir) >> 1;
              SWAP(arr[mid],arr[l+1]);
              if(arr[l+1] > arr[ir])
              {
                  SWAP(arr[l+1],arr[ir]);
              }
              if(arr[l] > arr[ir])
              {
                  SWAP(arr[l],arr[ir]);
              }
              if(arr[l+1] > arr[l])
              {
                  SWAP(arr[l+1],arr[l]);
              }
              i=l+1;        // initialize pointers for partitioning
              j=ir;
              a=arr[l];        // partitioning element
              for(;;)        // innermost loop
              {
                  do i++; while(arr[i]<a);        // scan up to find element > a
                  do j--; while(arr[j]>a);        // scan down to find element < a
                  if( j < i ) break;                // pointers crossed, median is in between
                  SWAP(arr[i],arr[j]);
              }
              arr[l]=arr[j];                        // insert partitioning element
              arr[j]=a;
              if( j >= k ) ir=j-1;                // Keep partition that contains the median active
              if( j <= k ) l=i;
          }
      }
  }
  
  
  double fmedian_with_rejection(unsigned int n, double arr[],double sigma_cutoff, double *final_mad, int *final_n)
  {
      double median_value;
      int i,done;
      double deviate,mad;
  
      done = 0;
      while(! done)
      {
          /* compute the median (centroid) value */
          median_value = fmedian(n,arr);
  
          /* now figure out what the mean absolute deviation from this value is */
          mad = fmedian_absolute_deviation(n,arr,median_value);
          //if(flag) printf("mad = %f\n",mad);
  
          done = 1;
          /* reject all outliers */
          for(i=1;i<=n;++i)
          {
              /* reject positive and negative outliers */
              deviate = fabs(arr[i]-median_value);
              if(deviate > sigma_cutoff*mad)
              {
                  /* needs to go */
                  /* move value at the end of the array to this "reject" and then shorten the array */
                  //if(flag) printf("rejecting arr[%d] = %f (%f)\n",i,arr[i],deviate);
                  //arr[worst]+=10000;
                  if(i != n)
                  {
                      //temp=arr[worst];
                      arr[i] = arr[n];
                      //arr[n]=temp;
                  }
                  --n;
                  done = 0;
              }
          }
      }
  
      /* basically three return values */
      *final_mad = mad;
      *final_n = n;
      return median_value;
  }
  
  /* note: there must be 2*n elements in this array! */
  double fmedian_absolute_deviation(unsigned int n, double arr[], double median_value)
  {
      int i;
      for(i=1;i<=n;++i)
      {
          arr[i+n] = fabs(arr[i]-median_value);
      }
  
      return fmedian(n,arr+n);
  }
  
  
  
  
  /* this function keeps track of outliers by swapping them to the end of the array */
  /* counting starts at 0 and "points" is the number of points */
  double fmean_with_rejection(unsigned int starting_points, double arr[], double sigma_cutoff, double *final_rmsd, int *final_n)
  {
      int points,i;
      int rejection,worst;
      double temp,sum,avg,sumd,rmsd,deviate,worst_deviate;
  
      points=starting_points;
  
      rejection = 1;
      avg = rmsd = NAN;
      while ( rejection && points>starting_points/2.0 )
      {
          /* find the mean and rms deivation */
          sum = sumd = 0.0;
          for(i=0;i<points;++i)
          {
              sum+=arr[i];
          }
          avg=sum/points;
          worst=-1;
          worst_deviate=0.0;
          for(i=0;i<points;++i)
          {
              deviate=fabs(arr[i]-avg);
              if(deviate > worst_deviate)
              {
                  worst=i;
                  worst_deviate=deviate;
              }
              sumd+=deviate*deviate;
          }
          rmsd=sqrt(sumd/points);
  
          rejection=0;
          if(worst_deviate>sigma_cutoff*rmsd)
          {
              /* we have a reject! */
              rejection=1;
  
              /* move it to end of the array and forget about it */
              SWAP(arr[worst],arr[points]);
              --points;
          }
      }
  
      *final_rmsd = rmsd;
      *final_n = points;
      return avg;
  }
  
  
  }}// namespace simtbx::nanoBragg
  