Page MenuHomec4science
Files F65983893

outerpotential.c
No OneTemporary
Actions

File Metadata

Created
Fri, Jun 7, 12:38

outerpotential.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include "allvars.h"
#include "proto.h"
#ifdef OUTERPOTENTIAL
/**************************************
*
* * * * driving functions * * *
*
***************************************/
void init_outer_potential()
{
#ifdef NFW
init_outer_potential_nfw();
#endif
#ifdef PLUMMER
init_outer_potential_plummer();
#endif
#ifdef MIYAMOTONAGAI
init_outer_potential_miyamotonagai();
#endif
#ifdef PISOTHERM
init_outer_potential_pisotherm();
#endif
#ifdef CORIOLIS
init_outer_potential_coriolis();
#endif
}
void outer_forces()
{
#ifdef NFW
outer_forces_nfw();
#endif
#ifdef PLUMMER
outer_forces_plummer();
#endif
#ifdef PISOTHERM
outer_forces_pisotherm();
#endif
#ifdef MIYAMOTONAGAI
outer_forces_miyamotonagai();
#endif
#ifdef CORIOLIS
outer_forces_coriolis();
#endif
}
void outer_potential()
{
#ifdef NFW
outer_potential_nfw();
#endif
#ifdef PLUMMER
outer_potential_plummer();
#endif
#ifdef MIYAMOTONAGAI
outer_potential_miyamotonagai();
#endif
#ifdef PISOTHERM
outer_potential_pisotherm();
#endif
#ifdef CORIOLIS
outer_potential_coriolis();
#endif
}
#ifdef NFW
/****************************
*
* * * * NFW * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_nfw()
{
double rhoc;
All.HaloMass = All.HaloMass / All.UnitMass_in_g;
All.HaloMass = All.HaloMass*SOLAR_MASS; /* from Msolar in user units */
All.NFWPotentialCte = (1.-All.GasMassFraction) * All.G*All.HaloMass/(log(1+All.HaloConcentration) - All.HaloConcentration/(1+All.HaloConcentration));
/* here, convert HubbleParam in 1/(unit_time) */
rhoc = pow(All.HubbleParam*HUBBLE*All.UnitTime_in_s,2)*3/(8*PI*All.G);
All.Rs = pow(( All.HaloMass/(4/3.*PI*200*rhoc) ),(1./3.)) / All.HaloConcentration;
printf("\nRs (internal units) = %g\n", All.Rs);
}
/*! compute forces du to the outer potential
*
*/
void outer_forces_nfw()
{
int i;
FLOAT accx, accy, accz;
FLOAT r,dPhi;
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
dPhi = - All.NFWPotentialCte * ( log(1+r/All.Rs) - r/(All.Rs+r) ) / (r*r);
accx = dPhi * P[i].Pos[0]/r;
accy = dPhi * P[i].Pos[1]/r;
accz = dPhi * P[i].Pos[2]/r;
/*
if (P[i].ID==8)
{
printf("------> x=%g y=%g z=%g ax=%g ay=%g az=%g r=%g\n",P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],accx,accy,accz,r);
printf("------> x=%g y=%g z=%g ax=%g ay=%g az=%g\n",P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],P[i].GravAccel[0],P[i].GravAccel[1],P[i].GravAccel[2]);
}
*/
P[i].GravAccel[0] = P[i].GravAccel[0] + accx;
P[i].GravAccel[1] = P[i].GravAccel[1] + accy;
P[i].GravAccel[2] = P[i].GravAccel[2] + accz;
}
}
}
/*! compute outer potential
*
*/
void outer_potential_nfw()
{
int i;
FLOAT r,Phi;
for(i = 0; i < NumPart; i++)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
Phi = - All.NFWPotentialCte * (log(1+r/All.Rs) / r);
/*
if (P[i].ID==8)
{
printf("------> x=%g y=%g z=%g Phi=%g\n",P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Phi);
printf("------> x=%g y=%g z=%g Phi=%g\n",P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],P[i].Potential);
}
*/
/* !!! do not forget do multiply by 2 (Springel choice ?) !!! */
P[i].Potential = P[i].Potential + 2*Phi;
}
}
#endif
#ifdef PLUMMER
/****************************
*
* * * * PLUMMER * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_plummer()
{
All.PlummerPotentialCte = All.G*All.PlummerMass /All.UnitMass_in_g *SOLAR_MASS; /* from Msolar in user units */
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_plummer()
{
int i;
FLOAT accx, accy, accz;
FLOAT r,dPhi;
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
dPhi = - All.PlummerPotentialCte / pow(r*r+All.PlummerSoftenning*All.PlummerSoftenning,1.5)*r;
accx = dPhi * P[i].Pos[0]/r;
accy = dPhi * P[i].Pos[1]/r;
accz = dPhi * P[i].Pos[2]/r;
P[i].GravAccel[0] = P[i].GravAccel[0] + accx;
P[i].GravAccel[1] = P[i].GravAccel[1] + accy;
P[i].GravAccel[2] = P[i].GravAccel[2] + accz;
}
}
}
/*! compute outer potential
*
*/
void outer_potential_plummer()
{
int i;
FLOAT r,Phi;
for(i = 0; i < NumPart; i++)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
Phi = - All.PlummerPotentialCte / pow(r*r+All.PlummerSoftenning*All.PlummerSoftenning,0.5);
/* !!! do not forget do multiply by 2 (Springel choice ?) !!! */
P[i].Potential = P[i].Potential + 2*Phi;
}
}
#endif
#ifdef MIYAMOTONAGAI
/****************************
*
* * * * MIYAMOTONAGAI * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_miyamotonagai()
{
All.MiyamotoNagaiPotentialCte = All.G*All.MiyamotoNagaiMass /All.UnitMass_in_g *SOLAR_MASS; /* from Msolar in user units */
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_miyamotonagai()
{
int i;
FLOAT accx, accy, accz;
FLOAT r,dPhi;
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
dPhi = - All.MiyamotoNagaiPotentialCte / pow(r*r+All.MiyamotoNagaiHr*All.MiyamotoNagaiHr,1.5)*r;
accx = dPhi * P[i].Pos[0]/r;
accy = dPhi * P[i].Pos[1]/r;
accz = dPhi * P[i].Pos[2]/r;
P[i].GravAccel[0] = P[i].GravAccel[0] + accx;
P[i].GravAccel[1] = P[i].GravAccel[1] + accy;
P[i].GravAccel[2] = P[i].GravAccel[2] + accz;
}
}
}
/*! compute outer potential
*
*/
void outer_potential_miyamotonagai()
{
int i;
FLOAT r,Phi;
for(i = 0; i < NumPart; i++)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
Phi = - All.MiyamotoNagaiPotentialCte / pow(r*r+All.MiyamotoNagaiHr*All.MiyamotoNagaiHr,0.5);
/* !!! do not forget do multiply by 2 (Springel choice ?) !!! */
P[i].Potential = P[i].Potential + 2*Phi;
}
}
#endif
#ifdef PISOTHERM
/**********************************************
*
* * * * PSEUDO ISOTHERM POTENTIAL * * *
*
***********************************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_pisotherm()
{
All.PisothermPotentialCte = 4*PI*All.Rho0*pow(All.Rc,3) *All.G *(1.-All.GasMassFraction);
All.Potentialw = gsl_integration_workspace_alloc (1000);
All.PotentialF.function = &potential_f;
/* potential normalisation */
/*All.PotentialInf = All.PisothermPotentialCte/All.Rc; viriel Pfen. */
All.PotentialInf = 3.3263051716357208; /* in order to have 2Uint = Epot */
}
/*! potential function
*
*/
double potential_f (double r, void * params)
{
double f = All.PisothermPotentialCte*(r/All.Rc-atan(r/All.Rc))/(r*r);
return f;
}
/*! potential function
*
*/
double get_potential(double r)
{
double result, error;
gsl_integration_qags (&All.PotentialF, 0, r, 0, 1e-5, 1000, All.Potentialw, &result, &error);
return result;
}
/*! compute forces du to the outer potential
*
*/
void outer_forces_pisotherm()
{
int i;
FLOAT accx, accy, accz;
FLOAT r,dPhi;
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
dPhi = - All.PisothermPotentialCte * ( r/All.Rc - atan(r/All.Rc) ) / (r*r);
accx = dPhi * P[i].Pos[0]/r;
accy = dPhi * P[i].Pos[1]/r;
accz = dPhi * P[i].Pos[2]/r;
P[i].GravAccel[0] = P[i].GravAccel[0] + accx;
P[i].GravAccel[1] = P[i].GravAccel[1] + accy;
P[i].GravAccel[2] = P[i].GravAccel[2] + accz;
}
}
}
/*! compute outer potential
*
*/
void outer_potential_pisotherm()
{
int i;
FLOAT Phi;
FLOAT r;
for(i = 0; i < NumPart; i++)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
Phi = get_potential(r)-All.PotentialInf;
/* !!! do not forget do multiply by 2 (Springel choice ?) !!! */
P[i].Potential = P[i].Potential + 2*Phi;
}
}
#endif
#ifdef CORIOLIS
/****************************
*
* * * * CORIOLIS * * *
*
*****************************/
/*! \file phase.c
* \brief Compute forces and potential from the outer potential
*
*/
/*! init outer potential constantes
*
*/
void init_outer_potential_coriolis()
{
set_outer_potential_coriolis();
}
void set_outer_potential_coriolis()
{
double amp;
double CoriolisTime = 1000;
amp = All.Time/CoriolisTime;
if (amp>1)
amp = 1;
All.CoriolisOmegaX = All.CoriolisOmegaX0*amp;
All.CoriolisOmegaY = All.CoriolisOmegaY0*amp;
All.CoriolisOmegaZ = All.CoriolisOmegaZ0*amp;
printf("coriolis CoriolisOmegaZ = %g\n",All.CoriolisOmegaZ);
}
/*! compute forces due to the outer potential
*
*/
void outer_forces_coriolis()
{
int i;
FLOAT accx, accy, accz;
for(i = 0; i < NumPart; i++)
{
if(P[i].Ti_endstep == All.Ti_Current)
{
accx = pow(All.CoriolisOmegaZ,2)*P[i].Pos[0] - All.CoriolisOmegaZ*P[i].Vel[1];
accy = pow(All.CoriolisOmegaZ,2)*P[i].Pos[1] + All.CoriolisOmegaZ*P[i].Vel[0];
accz = 0;
P[i].GravAccel[0] = P[i].GravAccel[0] + accx;
P[i].GravAccel[1] = P[i].GravAccel[1] + accy;
P[i].GravAccel[2] = P[i].GravAccel[2] + accz;
}
}
}
/*! compute outer potential
*
*/
void outer_potential_coriolis()
{
int i;
FLOAT r,Phi;
for(i = 0; i < NumPart; i++)
{
r = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
Phi = -0.5*pow(All.CoriolisOmegaZ*r,2);
/* !!! do not forget do multiply by 2 (Springel choice ?) !!! */
P[i].Potential = P[i].Potential + 2*Phi;
}
}
#endif
#endif

Event Timeline

c4science · Help