API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Nov 22nd, 2017
122
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 10.70 KB | None | 0 0
raw download clone embed print report
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <string.h>
  5. #define PI (3.1415926535897932384626433832795028841971693994)
  6. #include "init.h"
  7. #include "dens.h"
  8. #include "fields.h"
  9. #include "sources.h"
  10. #include "init.h"
  11.  
  12. struct fields *LOCAL_MESONS;
  13. struct sources *LOCAL_SRCS;
  14. struct densities *LOCAL_DENS;
  15. struct init_data COMMON_DATA;
  16.  
  17. //Need to add update Fields and Update Sigma then it's OVER!!!1!
  18.  
  19. double UpdateDensities(struct fields *Mesons, struct sources *srcs, struct densities *Dens, double mu_p, double mu_n, double rmax, int nmax);
  20. double UpdateProtonDensities();
  21. double UpdateNeutronDensities();
  22. double CombineOldAndNew();
  23. double solveNewton(double (*f)(double), double c, double x_init);
  24. double fPrime(double x, double (*f)(double));
  25. double CalcProtonNumber(double mu_p);
  26. double CalcNeutronNumber(double mu_n);
  27. void GetRho(double *vecRhopn, double mup, double *vecSigma, double *vecRho0, double *vecOmega, double *vecA0, int nMax, int p_or_n);
  28. double numIntegrateVec(double *vF, double h, int down, int up, int nMax);
  29. double GetEFermi(double mu, double eA, double grho, double gomega, int p_or_n);
  30. double GetPFermi(double EFermi, double mass);
  31. double ScalarDensity2(double m_eff, double pf);
  32. void Get_gRhoScalar(double *r_grho_S, double *sigma, double *rho_p, double *rho_n, double rmax, int nmax);
  33. double SinhXoverX(double x);
  34. void UpdateFields(struct fields *Mesons, struct sources *Srcs, struct densities *Dens, double rmax, int nmax);
  35. void UpdateSigma(struct fields *Mesons, struct densities *Dens, double rmax, int nmax);
  36. void SolveKG_WithSource(double *vecPhi, double *vecJ, double **matGreen, double m, double rMax, int nMax);
  37. void SolvePoisson(double *vecPhi, double *vecJ, double **matGreen, double rMax, int nMax);
  38. void MakeGreenE(double **matGreen, double rMax, int nMax);
  39.  
  40.  
  41.  
  42.  
  43.  
  44. //Update Densities
  45. double UpdateDensities(struct fields *Mesons, struct sources *SRCS, struct densities *Dens, double mu_p, double mu_n, double rmax, int nmax){
  46.   *LOCAL_MESONS = *Mesons;
  47.   *LOCAL_SRCS = *SRCS;
  48.   //Let LOCAL_DENS point to Dens.
  49.   *LOCAL_DENS = *Dens;
  50.  
  51.   UpdateProtonDensities(Dens->rho_p, mu_p);
  52.   //Pass by reference so both can be updated
  53.   double f_p;
  54.   f_p = CombineOldAndNew(Dens->rho_p_prev, Dens->rho_p, nmax);
  55.   // combine the new solution and the previous one
  56.   //to help convergence
  57.   // f_p is the relative difference between the new density
  58.   // and the previous one
  59.   UpdateNeutronDensities(Dens->rho_n, mu_n);
  60.   double f_n;
  61.   // pass by reference so that both can be updated
  62.   f_n = CombineOldAndNew(Dens->rho_n_prev, Dens->rho_n, nmax);
  63.   //Let LOCAL_MESONS,SRCS,DENS point to a NULL pointer so that nothing else can mess it up
  64.   return (f_p + f_n)/2;
  65. }
  66.  
  67.  // passed by
  68. double UpdateProtonDensities(double mu_p){
  69.   double Z = COMMON_DATA.Z;
  70.   double mu_init = mu_p;
  71.   double mu_now = solveNewton(CalcProtonNumber, Z, mu_init);
  72.   mu_p = mu_now;
  73. }
  74.  
  75. double UpdateNeutronDensities(double mu_n){
  76.   double N = COMMON_DATA.N;
  77.   double mu_init = mu_n;
  78.   double mu_now = solveNewton(CalcNeutronNumber, N, mu_init);
  79.   mu_n = mu_now;
  80. }
  81.  
  82. double CalcProtonNumber(double mu_p){
  83.   //double rmax = COMMON_DATA.rmax;
  84.   //int nmax = COMMON_DATA.nmax;
  85.  
  86.   //Get local Mesons
  87.   //Get rho_p from local dens
  88.   double *integp1;
  89.   integp1 = (double *) malloc((COMMON_DATA.nmax) * sizeof(double));
  90.  
  91.   GetRho(LOCAL_DENS->rho_p, mu_p, LOCAL_MESONS->sigma, LOCAL_MESONS->rho0, LOCAL_MESONS->omega, LOCAL_MESONS->A0, COMMON_DATA.nmax, 1);
  92.  
  93.   for(int l = 0; l < COMMON_DATA.nmax; l++){
  94.     double dR = COMMON_DATA.rmax/COMMON_DATA.nmax;
  95.     integp1[l] = 4*PI*pow(l*dR,2)*LOCAL_DENS->rho_p[l];
  96.   }
  97.   double Solp1 = numIntegrateVec(integp1, COMMON_DATA.rmax/COMMON_DATA.nmax, 0, COMMON_DATA.nmax, 100);
  98.   return Solp1;
  99. }
  100.  
  101. double CalcNeutronNumber(double mu_n){
  102.   //double rmax = COMMON_DATA.rmax;
  103.   //int nmax = COMMON_DATA.nmax;
  104.  
  105.   //Get local Mesons
  106.   //Get rho_p from local dens
  107.   double *integn1;
  108.   integn1 = (double *) malloc((COMMON_DATA.nmax) * sizeof(double));
  109.  
  110.   GetRho(LOCAL_DENS->rho_n, mu_n, LOCAL_MESONS->sigma, LOCAL_MESONS->rho0, LOCAL_MESONS->omega, LOCAL_MESONS->A0, COMMON_DATA.nmax, 0);
  111.  
  112.   for(int l = 0; l < COMMON_DATA.nmax; l++){
  113.     double dR = COMMON_DATA.rmax/COMMON_DATA.nmax;
  114.     integn1[l] = 4*PI*pow(l*dR,2)*LOCAL_DENS->rho_n[l];
  115.   }
  116.   double Soln1 = numIntegrateVec(integn1, COMMON_DATA.rmax/COMMON_DATA.nmax, 0, COMMON_DATA.nmax, COMMON_DATA.nmax);
  117.   return Soln1;
  118. }
  119.  
  120.  
  121.  
  122.  
  123. double CombineOldAndNew(double *rho_p_prev, double *rho_p, int nmax){
  124.   // CLARIFICTION: Add double new_frac and old_frac in your initial data
  125.   // structure. Set these in FillData. old_frac = 1.0 - new_frac.
  126.   // new_frac = 0.6 is a good place to start
  127.   //Get new_frac from COMMON_DATA.
  128.   //Get old_frac from COMMON_DATA.
  129.  
  130.   double err;
  131.   for(int k = 0; k< nmax; k++){
  132.     rho_p[k] = COMMON_DATA.new_frac*rho_p[k] + COMMON_DATA.old_frac*rho_p_prev[k];
  133.     //Caculate the relative error
  134.     err += sqrt( pow((rho_p[k] - rho_p_prev[k]),2.0)/pow((rho_p[k]),2.0) );
  135.   }
  136.   memcpy(rho_p_prev, rho_p, COMMON_DATA.nmax*sizeof(double));
  137.   return err;
  138. }
  139.  
  140. void Get_gRhoScalar(double *r_grho_S, double *sigma, double *rho_p, double *rho_n, double rmax, int nmax){
  141.   double dR = rmax/nmax;
  142.  
  143.   for(int k=0; k<=nmax; k++){
  144.     double r = k*dR;
  145.     double scalDen = ScalarDensity(COMMON_DATA.m_proton - COMMON_DATA.g_sigma*sigma[k], pow(rho_p[k]*3.0*PI*PI, (1.0/3.0)));
  146.     r_grho_S[k] = scalDen*COMMON_DATA.g_sigma*r;
  147.   }
  148. }
  149.  
  150. double ScalarDensity2(double m_eff, double pf){
  151.   double ef = sqrt(m_eff*m_eff + pf*pf);
  152.   double rho_s = m_eff*( pf*ef - m_eff*m_eff*asinh(pf/m_eff) )/2.0/PI/PI;
  153.   return rho_s;
  154. }
  155.  
  156. double SinhXoverX(double x){
  157.   if(x < 1e-8){
  158.     return 1.0 + pow(x,2.0)/6.0 + pow(x,4.0)/120.0;
  159.   }else return sinh(x)/x;
  160. }
  161.  
  162. void UpdateFields(struct fields *Mesons, struct sources *Srcs, struct densities *Dens, double rmax, int nmax){
  163.   UpdateSigma(Mesons, Dens, rmax, nmax);
  164.   //Use Solve KG for omega and rho0
  165.   SolveKG_WithSource(Mesons->omega, Srcs->r_grho_B, Mesons->GreenOmega, COMMON_DATA.m_omega, rmax, nmax);
  166.   SolveKG_WithSource(Mesons->rho0, Srcs->r_grho_3, Mesons->GreenRho0, COMMON_DATA.m_rho, rmax, nmax);
  167.   SolvePoisson(Mesons->A0, Srcs->r_grho_p,Mesons->GreenA0, rmax, nmax);
  168.  
  169.   //Use SolvePossion for A0
  170.  
  171.   return;
  172. }
  173.  
  174. void UpdateSigma(struct fields *Mesons, struct densities *Dens, double rmax, int nmax){
  175.   double *sigma_prev = Mesons->sigma_prev = makeVector(com_dat.nmax + 1);
  176.   double *r_grho_S = Mesons->r_grho_S = makeVector(com_dat.nmax + 1);
  177.  
  178.   double *sigma;
  179.   sigma = Mesons->sigma;
  180.   double **GreenSigma = Mesons->GreenSigma;
  181.  
  182.   double *rho_p = Dens->rho_p;
  183.   double *rho_n = Dens->rho_n;
  184.  
  185.   double mass = COMMON_DATA.m_sigma;
  186.   //Get g_sigma from common_data
  187.  
  188.   int iter_max = 100;
  189.   double tol = 1e-10;
  190.   double err;
  191.   for(int i = 0; i<iter_max;){
  192.     memcpy(sigma_prev, sigma, COMMON_DATA.nmax*sizeof(double));
  193.     Get_gRhoScalar(r_grho_S, sigma, rho_p, rho_n, rmax, nmax);
  194.     SolveKG_WithSource(sigma, r_grho_S, GreenSigma, mass, rmax, nmax);
  195.  
  196.     //rho_p[k] = COMMON_DATA.new_frac*rho_p[k] + COMMON_DATA.old_frac*rho_p_prev[k];
  197.     //Caculate the relative error
  198.     err += sqrt( pow((sigma[i] - sigma_prev[i]),2.0)/pow((sigma[i]),2.0) );
  199.     if(err>tol){
  200.       break;
  201.     }
  202.   }
  203.   // Free r_g_rho_S and  sigma_prev;
  204.   return;
  205. }
  206.  
  207.  
  208.  
  209. //From previous assignment
  210. void GetRho(double *vecRhopn, double mup, double *vecSigma, double *vecRho0, double *vecOmega, double *vecA0, int nMax, int p_or_n){
  211.  
  212.   //Proton case
  213.   if(p_or_n == 1){
  214.     for (int k = 0; k < nMax; k++) {
  215.       //Make sure to use effective mass as given in assignment AND account for hbarc
  216.       double effMass = com_dat.m_proton/com_dat.hbarc - com_dat.g_sigma*vecSigma[k];
  217.       double eFerm = GetEFermi(mup, com_dat.e_EM*vecA0[k], com_dat.g_rho*vecRho0[k], com_dat.g_omega*vecOmega[k], p_or_n);
  218.       double pFerm = GetPFermi(eFerm, effMass);
  219.       //Populating the vector in given in the first argument.
  220.       vecRhopn[k] = pow(pFerm, 3)/(3.*pow(PI, 2));
  221.     }
  222.  
  223.   }
  224.   //Neutron case
  225.   else{
  226.     for (int k = 0; k < nMax; k++) {
  227.       double effMass = com_dat.m_neutron/com_dat.hbarc - com_dat.g_sigma*vecSigma[k];
  228.       double eFerm = GetEFermi(mup, com_dat.e_EM*vecA0[k], com_dat.g_rho*vecRho0[k], com_dat.g_omega*vecOmega[k], p_or_n);
  229.       double pFerm = GetPFermi(eFerm, effMass);
  230.       vecRhopn[k] = pow(pFerm, 3)/(3.*pow(PI, 2));
  231.     }
  232.   }
  233.  
  234. }
  235.  
  236.  
  237. double fPrime(double x, double (*f)(double)){
  238.   double h = 0.000000000000001;
  239.   return (f(x-3.0*h) - 27.0*f(x-h) +27.0*f(x+h) - f(x-3.0*h))/48.0*h;
  240. }
  241.  
  242. double solveNewton(double (*f)(double), double c, double x_init) {
  243.   //set the tolerance
  244.   double tol = 0.000000000000001;
  245.   //  set the maximum count (how many iterations)
  246.   int maxCount = 100000;
  247.   //-----dasdsalkjdlsajd
  248.   int n = 10;
  249.   //initializing x_old to be zero to start
  250.   double x_o = 0;
  251.   //Setting x_old to x_initial given in the function as an argument.
  252.   x_o = x_init;
  253.   //set x_new to be x_old with the n*tolerance added
  254.   double x_n = x_o + n*tol;
  255.   //initialize my counter to be zero to begin
  256.   int count = 0;
  257.   while (fabs(x_n - x_o)>tol){
  258.  
  259.     //Set x_old to x_new
  260.     x_o = x_n;
  261.     //Now make x_new = x_old -
  262.     x_n = x_o - ((*f)(x_o)-c)/fPrime(x_o, f);
  263.  
  264.     //count + 1
  265.     count += 1;
  266.  
  267.     //Testing
  268.     ///////////////////////////////////////////////////////////////////////////////////////
  269.     //printf("Iteration: %d, value is %.12f.\n", count, x_n);
  270.     //////////////////////////////////////////////////////////////////////////////////////
  271.  
  272.     //Break the loop if the count exceeds the maximum count set.
  273.     if (count > maxCount){
  274.       break;
  275.     }
  276.  
  277.     }
  278.     //The whole function now returns the latest iteration of x_new
  279.     return x_n;
  280. }
  281.  
  282. //Getting the fermi energy
  283. double GetEFermi(double mu, double eA, double grho, double gomega, int p_or_n){
  284.   //proton case
  285.   if(p_or_n == 1){
  286.     return mu - eA - grho - gomega;
  287.   }
  288.   //neutron case
  289.   else{
  290.     return mu + grho - gomega;
  291.   }
  292. }
  293. //Fermi momentum
  294. double GetPFermi(double EFermi, double mass){
  295.   return sqrt(pow(EFermi, 2) - pow(mass, 2));
  296. }
  297.  
  298. void SolveKG_WithSource(double *vecPhi, double *vecJ, double **matGreen, double m, double rMax, int nMax){
  299.   double dR = rMax/nMax;
  300.   //Variable we're going to use.
  301.   double dX = dR*m;
  302.  
  303.   double f_1, f_2;
  304.  
  305.   double *integ_vec;
  306.   integ_vec = (double *) malloc((nMax) * sizeof(double));
  307.  
  308.   for(int k=0; k<nMax; k++){
  309.     for(int l=0; l<nMax; l++){
  310.       integ_vec[l] = matGreen[k][l]*vecJ[l];
  311.     }
  312.  
  313.     f_1 = numIntegrateVec(integ_vec, dX, 0, k, nMax);
  314.     f_2 = numIntegrateVec(integ_vec, dX, k, nMax, nMax);
  315.     //Populating the solution vector
  316.     vecPhi[k] = (f_1 + f_2)/m;
  317.   }
  318.   return;
  319. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!