deprecated / bimsim

  #include "nearfield.h"

  #include <time.h>
  #include <math.h>
  
  #ifdef _WIN32
  #define isnan(x) _isnan(x)

  #define isinf(x) (!_finite(x))
  #endif
  
  int bessjyv_sph(int v, double z, double &vm, double* cjv,

      double* cyv, double* cjvp, double* cyvp);

  
  nearfieldStruct::nearfieldStruct()
  {
      scalarSim = true;
  	planeWave = false;

  	lut_us = true;
  	lut_uf = false;

  
  	nWaves = 0;
  }
  
  void nearfieldStruct::init()
  {
  	//set the field parameters
  	U.scalarField = scalarSim;
  	Uf.scalarField = scalarSim;
  
  	//initialize dynamic memory
  	U.init_gpu();
  	Uf.init_gpu();
  }
  
  void nearfieldStruct::destroy()
  {
  	U.kill_gpu();
  	Uf.kill_gpu();
  }
  
  void nearfieldStruct::setPos(bsPoint pMin, bsPoint pMax, bsVector normal)
  {
  	pos = rts::rtsQuad<ptype, 3>(pMin, pMax, normal);
  }
  
  void nearfieldStruct::setRes(int x_res, int y_res)
  {
  	U.R[0] = Uf.R[0] = x_res;
  	U.R[1] = Uf.R[1] = y_res;
  }
  
  std::string nearfieldStruct::toStr()
  {
  	std::stringstream ss;
  
  	ss<<"------Field Parameters-------"<<std::endl;
  	ss<<"Wavelength: "<<lambda<<"um"<<std::endl;
  	ss<<"K Vector (r, theta, phi): "<<k.cart2sph()<<std::endl;
  	ss<<"Condenser NA: "<<condenser[0]<<" to "<<condenser[1]<<std::endl;
  	ss<<"Focal Point: "<<focus[0]<<", "<<focus[1]<<", "<<focus[2]<<std::endl;
  	ss<<"Field Slice: "<<std::endl;

  	if(lut_us)
          ss<<"LUT Parameters --- min: "<<d_min<<"   max: "<<d_max<<std::endl;

  	ss<<pos<<std::endl;
  
  	ss<<std::endl<<"---------Materials-----------"<<std::endl;
  	ss<<"Number of Materials: "<<mVector.size()<<std::endl;
  	ss<<"Refractive Indices at lambda = "<<lambda<<"um"<<std::endl;
  	//output each material
  	for(unsigned int m=0; m<mVector.size(); m++)
  		ss<<" "<<m<<": "<<mVector[m](lambda)<<std::endl;
  
  	ss<<"---------Spheres-------------"<<std::endl;
  	ss<<"Number of Spheres: "<<sVector.size()<<std::endl;
  	//output each sphere
  	for(unsigned int s=0; s<sVector.size(); s++)
  		ss<<sVector[s].toStr()<<std::endl;
  

      ss<<"---------Timings-------------"<<std::endl;
      ss<<"Uf = "<<t_Uf<<"ms"<<std::endl;
      ss<<"Us = "<<t_Us<<"ms"<<std::endl;
  

  	return ss.str();
  }
  
  //generate monte-carlo waves
  void nearfieldStruct::calcWaves()
  {
      inWaves.resize(nWaves);
  
      //re-seed the random number generator

      //srand(time(NULL));
  	srand(NULL);

  
      //calculate the monte-carlo samples
      mcSampleNA(&inWaves[0], nWaves, k, condenser[0], condenser[1]);
  
  
  }
  
  
  void nearfieldStruct::calcSpheres()
  {
      //calculate all of the constants necessary to evaluate the scattered field
  	//estimate the order required to represent the scattered field for each sphere
  

  	//for each sphere

  	for(unsigned int i=0; i<sVector.size(); i++)

  	{
  		//a = sVector[i].a;
  
  		//calculate the required order
  		sVector[i].calcNl(lambda);

  
  		//set the refractive index for the sphere
  		int imat = sVector[i].iMaterial;
          rts::rtsComplex<ptype> n = mVector[imat](lambda);

  
  		//calculate the scattering coefficients
  		sVector[i].calcCoeff(lambda, n);
  
  		//save the refractive index
  		sVector[i].n = n;
  

  		//if the LUT is used, calculate Usp(theta, r)
  		if(lut_us)
  		{
  			sVector[i].calcUp(lambda, n, pos, max(U.R[0], U.R[1]));
          }
  
  

  	}
  
  }
  

  void nearfieldStruct::calcUs()
  {
  
  
      if(lut_us)
          scalarUpLut();
      else
          scalarUs();
  }
  
  void nearfieldStruct::calcUf()
  {
  	if(lut_uf)
  		scalarUfLut();
  	else
  		scalarUf();
  }
  

  void nearfieldStruct::Simulate()
  {

      //initialize timings
      t_Uf = 0;
      t_Us = 0;
  

  	//compute a set of plane waves for Monte-Carlo simulation
  	calcWaves();
  
      //the near field has to be simulated no matter what the output rtsPoint is

      calcUf();

      calcSpheres();

      calcUs();

      sumUf();

  	if(verbose)
          U.Mag().toImage("testU.bmp");

  }
  
  void nearfieldStruct::calcBesselLut(ptype* j, ptype d_min, ptype d_max, int dR)
  {
      /*Compute the look-up-table for spherical bessel functions used for the incident field
          j    =   (Nl + 1) x aR array of values
          aR      =   resolution of j
      */
  
  	//compute the wavenumber
  	ptype k = 2 * PI / lambda;
  	unsigned int Nl = m;
  
      //allocate space for the Bessel functions of the first and second kind (and derivatives -- which will be ignored)
      int bytes = sizeof(double) * (Nl + 1);
      double* cjv_kd = (double*)malloc(bytes);
      double* cyv_kd = (double*)malloc(bytes);
      double* cjvp_kd = (double*)malloc(bytes);
      double* cyvp_kd = (double*)malloc(bytes);
  
      //compute the bessel functions using the CPU-based algorithm
      double vm;
  
      //for each sample along r
      ptype dr = (d_max - d_min) / (dR - 1);
      ptype d;
      ptype jv;
      for(int id = 0; id < dR; id++)
      {
          d = id * dr + d_min;
          double kd = k*d;
          bessjyv_sph(Nl, kd, vm, cjv_kd, cyv_kd, cjvp_kd, cyvp_kd);
  
          //copy the double data to the bsComplex array

          for(unsigned int l=0; l<=Nl; l++)

  		{
              jv = cjv_kd[l];
  			if(isnan(jv) || isinf(jv))
  			{
                  if(kd == 0 && l == 0)
                      jv = 1;
                  else
                      jv = 0;
              }
              j[id * (Nl+1) + l] = jv;
  		}
      }
  

      if(verbose)

      {

          ofstream outfile("uf_besselout.txt");
          for(int ir = 0; ir < dR; ir++)

          {

              outfile<<ir*dr + d_min<<endl;

              for(unsigned int l = 0; l<=Nl; l++)

              {
                  outfile<<j[ir * (Nl+1) + l]<<" -- ";
              }
              outfile<<endl;

          }

          outfile.close();

      }

  }