/*
Soud Kharusi 260620452
PHYS 557 Assignment 7
Due Nov 22 2017
*/

//To compile
//gcc "filename.c" -c -o -lm execfilename
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define PI (3.1415926535897932384626433832795028841971693994)

#include "init.h"
#include "dens.h"
#include "fields.h"
#include "sources.h"
#include "init.c"
#include "update.c"

struct init_data com_dat;
struct fields mesons;
struct densities dens;
struct sources srcs;

//struct fields *LOCAL_MESONS;
//struct sources *LOCAL_SRCS;
//struct densities *LOCAL_DENS;
//struct init_data COMMON_DATA;

//declaring my functions to be used in main. Full functions are below main.
//p_or_n is 1 for proton and 0 for neutron
void FillData();
//double GetEFermi(double mu, double eA, double grho, double gomega, int p_or_n);
//double GetPFermi(double EFermi, double mass);
//void GetRho(double *vecRhopn, double mup, double *vecSigma, double *vecRho0, double *vecOmega, double *vecA0, int nMax, int p_or_n);


//The main part of the programme, where I will run and print as well as output the dat file.
int main() {
  //filling in all the initial data.
  FillData();
  //Allocate memory space for the density vectors and fill with initial guess.
  InitializeDensities(&dens);

  //Allocate memory space to meson (and A0) vectors and also make the Green function matrices
  InitializeFields(&mesons);

  //Allocate memory space for the source vectors and fill them with the initial guess
  InitializeSources(&srcs, &dens, &mesons);

  // Fill the meson and A0 vectors with the initial values according to the intial density values
  MakeInitFields(&mesons, &srcs, &dens);

  double mu_p = COMMON_DATA.m_proton; // Initial guess for the proton chemical potential
  double mu_n = COMMON_DATA.m_neutron; // Ditto for the neutrons
  double err = UpdateDensities(&mesons,&srcs,&dens, mu_p, mu_n, COMMON_DATA.rmax, COMMON_DATA.nmax);
  int iter_max = 100; // suitably large number. 100 will do
  for(int iter=0; iter<iter_max; iter++){
    UpdateSources(&srcs, &dens, &mesons, COMMON_DATA.rmax, COMMON_DATA.nmax);
    UpdateFields(&mesons, &srcs, &dens, COMMON_DATA.rmax, COMMON_DATA.nmax);
    err = UpdateDensities(&mesons,&srcs,&dens, mu_p, mu_n, COMMON_DATA.rmax, COMMON_DATA.nmax);
    double tol = 1e-10;
    if(err < tol){
      break;
    }
    //PrintOut(Mesons, Dens) // Print out densities and fields
    return 0;
  }
/*

  //Tests below
  double* integ = makeVector(com_dat.nmax + 1);
  double dr = com_dat.rmax/com_dat.nmax;
  for(int i = 0; i<com_dat.nmax+1; i++){
    integ[i] = 4.0*M_PI*pow(i*dr, 2)*dens.rho_n[i];
    //printf("%.12f\n",integ[i]);
  }
  double test = numIntegrateVec(integ, com_dat.rmax/com_dat.nmax , 0, com_dat.nmax, com_dat.nmax);

  printf("Test should give either Z=20 or N=22: %.12f \n", test);

*/

/*

  //Outputting the results to a .dat file.
  FILE *fp;
  // "w" means that we are going to write on this file
  fp = fopen("Kharusi_Result_HW6_labelled.dat", "w");
  double dr = com_dat.rmax/com_dat.nmax;
  for(int i = 0; i < com_dat.nmax; i++){
    fprintf(fp, "radius = %.2f \t sigma = %.12f \t A0 = %.12f \t omega = %.12f \t rho0= %.12f \n", i*dr, mesons.sigma[i]*com_dat.hbarc, mesons.A0[i]*com_dat.hbarc, mesons.omega[i]*com_dat.hbarc, mesons.rho0[i]*com_dat.hbarc);
  }
  fclose(fp);

  FILE *fp2;
  // "w" means that we are going to write on this file
  fp2 = fopen("Kharusi_Result_HW6.dat", "w");
  for(int i = 0; i < com_dat.nmax; i++){
    fprintf(fp2, "%.2f \t %.12f \t %.12f \t %.12f \t %.12f \n", i*dr, mesons.sigma[i]*com_dat.hbarc, mesons.A0[i]*com_dat.hbarc, mesons.omega[i]*com_dat.hbarc, mesons.rho0[i]*com_dat.hbarc);
  }
  fclose(fp2);
*/

return 0;
}



void FillData(){
  //Example for this assignment
  com_dat.nmax = 200;
  com_dat.rmax = 12.0;
  com_dat.A = 42.;
  com_dat.Z = 20.;
  com_dat.N = 22.;

  com_dat.hbarc = 197.3;
  com_dat.m_rho = 770./com_dat.hbarc;
  com_dat.m_omega = 783./com_dat.hbarc;
  com_dat.m_sigma = 550./com_dat.hbarc;
  com_dat.g_rho = 2.629;
  com_dat.e_EM = 0.3028221199;
  com_dat.g_omega = 13.80557973;
  com_dat.g_sigma = 11.07849711;
  com_dat.m_proton = 938.27/com_dat.hbarc;
  //com_dat.m_proton = 939.0/com_dat.hbarc;
  com_dat.m_neutron = 939.57/com_dat.hbarc;
  //com_dat.m_neutron = 939.0/com_dat.hbarc;
  COMMON_DATA.old_frac = 0.6;
  COMMON_DATA.new_frac = 1.0 - COMMON_DATA.old_frac;
}