#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <complex>
#include <random>

const double PI=4.0*atan(1.0);
const double HBARC=197.3269602;

using namespace std;

int main(int argc,char *argv[]){
	
	// Here is an example of complex arithmetic
  complex<double> ci(0.0,1.0);
  complex<double> z;
  double phi;
  printf("Enter phi in degrees: ");
  scanf("%lf",&phi);
  z=exp(ci*phi*PI/180.0);
  printf("phi=%g, e^iphi=(%g,%g)\n",phi,real(z),imag(z));
	
	// Here is an example of some loops
	int n;
	double sum=0,piguess;
	for(n=1;n<=1000;n++){
		sum=sum+1.0/pow(double(n),4);
	}
	piguess=pow(90.0*sum,0.25);
	printf("PI=%12.9f =? %12.9f\n",PI,piguess);
	
	// You can do the same with loops
	sum=0; n=1;
	do{
		sum+=1.0/pow(double(n),4);
		n+=1;
	} while(n<=1000);
	printf("PI=%12.9f =? %12.9f\n",PI,piguess);
	sum=0; n=1;
	while(n<=1000){
		sum+=1.0/pow(double(n),4);
		n+=1;
	}
	printf("PI=%12.9f =? %12.9f\n",PI,piguess);
	
	// Random practice
	int seed=time(NULL);
	mt19937 rng(seed);
	std::normal_distribution<double> gauss_dist(0,1);
	double x,x2bar=0.0;
	int nmax=10000000;
	for(n=0;n<nmax;n++){
		x=gauss_dist(rng);
		x2bar+=x*x;
	}
	x2bar=x2bar/double(nmax);
	printf("x2bar=%g =? 1.0\n",x2bar);
	
	
  return 0;
}