/*
 * Seminarski rad iz Fizike3 - Fizika materijala
 * prof. Horvat
 * autor U*nix porta: Dinko Korunic, 36355514, email: kreator@fly.srk.fer.hr
 * autor DOS porta: Domagoj Kristo, xxx, email: Domagoj.Kristo@fer.hr
 *
 * Zadatak: Izracunati faktor transmisije T da elektron tunelira kroz
 * potencijalnu barijeru (pomocu WBK aproksimacije) koja je zadana sa
 * V[x_]=[x/10^-14]*6MeV. Doticna vjerojatnost jednaka je (pri pocetnim
 * uvjetima da je Emax=6MeV i da je Epotrebna=E=3MeV)
 * T=exp[-2Sqrt[4*pi*me/h]Integrate[Sqrt[V[x]-E]],{x,z,10^-14}]. Dobiveni
 * pocetni uvjet z=NSolve[V[x]==E,x] a to je z=2^-1/3*10^-14. Za detaljniji
 * opis vidjeti prilozeti Mathematica3 dokument.
 *
 * Doticnu integraciju treba izvesti MonteCarlo metodom koja se bazira na
 * slijedecem: f-ja je y(x) i zatvorimo ju sa pravokutnikom sirine b-a i
 * visinom c, dakle povrsine Ar=(b-c)a. Onda generiramo slucajne tocke
 * _unutar_ doticnog pravokutnika (xr,yr). Onda izracunamo da li je nasa tocka
 * pogodak ili promasaj time sto za xr izracunamo f(xr) i usporedimo sa yr.
 * Izracunavamo mnogo takvih tocaka povecavajuci brojac pogodaka i brojac
 * pokusaja.. Trazena povrsina je (broj_pogodaka/broj_pokusaja)*Ar. Mozemo
 * nastaviti sa pokusajima dok se ne postigne zeljena preciznost.
 *
 * Please enter resolution: 99999999
 * Calculus values:
 * x1=7.937005e-15
 * x2=1.000000e-14
 * height=1.749371e+00
 * width=2.062995e-15
 * area A=3.608944e-15
 * Calculating...
 * trials 99999999, hits 62894496
 * Result for 99999999 iterations would be 2.269827e-15
 *
 * $Id: mcarl.c,v 1.2 1999/11/17 15:33:09 kreator Exp kreator $
 * */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>

/* integrabilna f-ja */
#define mydeffunction(x) sqrt(pow((x/1E-14),3.)*6.-3.)
/* donja granica integracije */
#define mydefa pow(2.,(-1./3.))*1E-14
/* gornja granica integracije */
#define mydefb 1E-14 
/* approx visina pravokutnika, moramo se osigurati da bude veca od max f(x) */
#define mydefc 1.001*myfunction(myb)

double myfunction(double x)
{
	return mydeffunction(x);
}

double myintegrate(long *myprecision)
{
	long trials,		/* broj pokusaja */
		hits;		/* broj pogodaka */
	double Ar,		/* approx povrsina */
		A,		/* trazena povrsina */
		mya=mydefa, 	/* donja */
		myb=mydefb, 	/* gornja */
		myc=mydefc,	/* approx visina */
		myrndx,		/* random x */
		myrndy; 	/* random y*/

	/* approx povrsina */
	Ar=myc*(myb-mya);
	/* random seed */
	srandom(getppid()+getpid()+time(NULL));

	/* be more verbose */	
	fprintf(stdout, "Calculus values:\nx1=%e\nx2=%e\nheight=%e\n"
			"width=%e\narea A=%e\n", mya, myb, myc, myb-mya, Ar);
	
	fprintf(stdout, "Calculating...\n");
	for (trials=hits=0; trials<*myprecision; ++trials)
	{
		myrndx=random()/(RAND_MAX+1.)*(myb-mya)+mya;
		/* if (myrndx>myb || myrndx<mya) exit(1); */
		myrndy=random()/(RAND_MAX+1.)*myc;
		if (myfunction(myrndx)>=myrndy)
			/* pogodak! */
			++hits;
	}

	fprintf(stdout, "trials %ld, hits %ld\n", trials, hits);
	return A=Ar*hits/trials;
}

int main(void)
{
	long myresolution;
	double myresult;

	fprintf(stdout, "Please enter resolution: ");
	fscanf(stdin, "%ld", &myresolution);
	myresult=myintegrate(&myresolution);
	fprintf(stdout, "Result for %ld iterations would be %e\n",
			myresolution, myresult);
	return EXIT_SUCCESS;
}