testing/icgen/random_verschillende_resoluties_N-GenIC/gadget3_64/get_Je.c
int star_density_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int j, n, numngb;
  int startnode, listindex = 0;
  double h, hinv, h2, mass_j, weight, hinv3;
  double wk, mass, density, lum;
  double dx, dy, dz, r, r2, u, a3inv;
  MyDouble *pos;

#ifdef PERIODIC
  double boxsize, boxhalf;

  boxsize = All.BoxSize;
  boxhalf = 0.5 * All.BoxSize;
#endif

  if(All.ComovingIntegrationOn)
    a3inv = 1.0 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1.0;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
      density = P[target].DensAroundStar;
      mass = P[target].Mass;
    }
  else
    {
      pos = StarDataGet[target].Pos;
      h = StarDataGet[target].Hsml;
      density = StarDataGet[target].Density;
      mass = StarDataGet[target].Mass;
    }

  h2 = h * h;
  hinv = 1.0 / h;
  hinv3 = hinv * hinv * hinv;


  lum = mass * All.IonizingLumPerSolarMass * (PROTONMASS / SOLAR_MASS) * All.UnitTime_in_s / All.HubbleParam;

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = StarDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  numngb = ngb_treefind_variable(pos, h, target, &startnode, mode, nexport, nsend_local);

	  if(numngb < 0)
	    return -1;

	  for(n = 0; n < numngb; n++)
	    {
	      j = Ngblist[n];

	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];
#ifdef PERIODIC
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;
	      r = sqrt(r2);

	      if(r2 < h2)
		{
		  u = r * hinv;

		  if(u < 0.5)
		    wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		  else
		    wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);

		}
	      else
		wk = 0;

	      mass_j = P[j].Mass;

	      weight = mass_j * wk / density;

	      SphP[j].Je += lum * weight / mass_j;
	    }
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = StarDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  return 0;
}
/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
This function represents the core of the SPH density computation. The
target particle may either be local, or reside in the communication
buffer.
This function represents the core of the SPH density computation . The target particle may either be local or reside in the communication buffer .
int star_density_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int j, n, numngb;
  int startnode, listindex = 0;
  double h, hinv, h2, mass_j, weight, hinv3;
  double wk, mass, density, lum;
  double dx, dy, dz, r, r2, u, a3inv;
  MyDouble *pos;
#ifdef PERIODIC
  double boxsize, boxhalf;
  boxsize = All.BoxSize;
  boxhalf = 0.5 * All.BoxSize;
#endif
  if(All.ComovingIntegrationOn)
    a3inv = 1.0 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1.0;
  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
      density = P[target].DensAroundStar;
      mass = P[target].Mass;
    }
  else
    {
      pos = StarDataGet[target].Pos;
      h = StarDataGet[target].Hsml;
      density = StarDataGet[target].Density;
      mass = StarDataGet[target].Mass;
    }
  h2 = h * h;
  hinv = 1.0 / h;
  hinv3 = hinv * hinv * hinv;
  lum = mass * All.IonizingLumPerSolarMass * (PROTONMASS / SOLAR_MASS) * All.UnitTime_in_s / All.HubbleParam;
  if(mode == 0)
    {
      startnode = All.MaxPart;	
    }
  else
    {
      startnode = StarDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	
    }
  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  numngb = ngb_treefind_variable(pos, h, target, &startnode, mode, nexport, nsend_local);
	  if(numngb < 0)
	    return -1;
	  for(n = 0; n < numngb; n++)
	    {
	      j = Ngblist[n];
	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];
#ifdef PERIODIC
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;
	      r = sqrt(r2);
	      if(r2 < h2)
		{
		  u = r * hinv;
		  if(u < 0.5)
		    wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		  else
		    wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
		}
	      else
		wk = 0;
	      mass_j = P[j].Mass;
	      weight = mass_j * wk / density;
	      SphP[j].Je += lum * weight / mass_j;
	    }
	}
      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = StarDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	
	    }
	}
    }
  return 0;
}
int star_density_evaluate ( int target , int mode , int * nexport , int * nsend_local ) { int j , n , numngb ; int startnode , listindex = 0 ; double h , hinv , h2 , mass_j , weight , hinv3 ; double wk , mass , density , lum ; double dx , dy , dz , r , r2 , u , a3inv ; MyDouble * pos ; #ifdef PERIODIC double boxsize , boxhalf ; boxsize = All . BoxSize ; boxhalf = 0.5 * All . BoxSize ; #endif if ( All . ComovingIntegrationOn ) a3inv = 1.0 / ( All . Time * All . Time * All . Time ) ; else a3inv = 1.0 ; if ( mode == 0 ) { pos = P [ target ] . Pos ; h = PPP [ target ] . Hsml ; density = P [ target ] . DensAroundStar ; mass = P [ target ] . Mass ; } else { pos = StarDataGet [ target ] . Pos ; h = StarDataGet [ target ] . Hsml ; density = StarDataGet [ target ] . Density ; mass = StarDataGet [ target ] . Mass ; } h2 = h * h ; hinv = 1.0 / h ; hinv3 = hinv * hinv * hinv ; lum = mass * All . IonizingLumPerSolarMass * ( PROTONMASS / SOLAR_MASS ) * All . UnitTime_in_s / All . HubbleParam ; if ( mode == 0 ) { startnode = All . MaxPart ; } else { startnode = StarDataGet [ target ] . NodeList [ 0 ] ; startnode = Nodes [ startnode ] . u . d . nextnode ; } while ( startnode >= 0 ) { while ( startnode >= 0 ) { numngb = ngb_treefind_variable ( pos , h , target , & startnode , mode , nexport , nsend_local ) ; if ( numngb < 0 ) return -1 ; for ( n = 0 ; n < numngb ; n ++ ) { j = Ngblist [ n ] ; dx = pos [ 0 ] - P [ j ] . Pos [ 0 ] ; dy = pos [ 1 ] - P [ j ] . Pos [ 1 ] ; dz = pos [ 2 ] - P [ j ] . Pos [ 2 ] ; #ifdef PERIODIC if ( dx > boxHalf_X ) dx -= boxSize_X ; if ( dx < - boxHalf_X ) dx += boxSize_X ; if ( dy > boxHalf_Y ) dy -= boxSize_Y ; if ( dy < - boxHalf_Y ) dy += boxSize_Y ; if ( dz > boxHalf_Z ) dz -= boxSize_Z ; if ( dz < - boxHalf_Z ) dz += boxSize_Z ; #endif r2 = dx * dx + dy * dy + dz * dz ; r = sqrt ( r2 ) ; if ( r2 < h2 ) { u = r * hinv ; if ( u < 0.5 ) wk = hinv3 * ( KERNEL_COEFF_1 + KERNEL_COEFF_2 * ( u - 1 ) * u * u ) ; else wk = hinv3 * KERNEL_COEFF_5 * ( 1.0 - u ) * ( 1.0 - u ) * ( 1.0 - u ) ; } else wk = 0 ; mass_j = P [ j ] . Mass ; weight = mass_j * wk / density ; SphP [ j ] . Je += lum * weight / mass_j ; } } if ( mode == 1 ) { listindex ++ ; if ( listindex < NODELISTLENGTH ) { startnode = StarDataGet [ target ] . NodeList [ listindex ] ; if ( startnode >= 0 ) startnode = Nodes [ startnode ] . u . d . nextnode ; } } } return 0 ; }
This function represents the core of the SPH density computation.
This function represents the core of the SPH density computation .
/* root node */ /* open it */ /* open it */