mc.cpp

用Monte Carlo方法

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

#include "Reaction.h"

#define STEP 10	//输出间隔

double stime=0.0;

//子函数
void print(FILE *fp);	//输出变量值
void go();				//进行一次抽样

void main()
{
	int i=0;
	FILE *fp;

	fp=fopen("out.txt","w");
	init();
	fprintf(fp,id_parameter);

	do
	{
		for(i=0;i<STEP;i++)
		{
			equilibrium();
			go();
		}
		print(fp);
	}while(stime<40);

	T=42;
	do
	{
		for(i=0;i<STEP;i++)
		{
			equilibrium();
			go();
		}
		print(fp);
	}while(stime<400);
	fclose(fp);
}

void print(FILE *fp)
{
	fprintf(fp,"%e\t",stime);
	for(int i=0;i<n_parameter;i++)fprintf(fp,"%e\t",*parameter[i]);
	fprintf(fp,"\n");
}

void go()
{
	double pa[n_channel+1];
	double pi[n_channel+1];
	int i,j;

	pa[0]=0.0;pi[0]=0.0;
	for(i=1;i<=n_channel;i++)
	{
		pa[i]=poss[i]();pi[i]=0;
	}


	for(i=1;i<=n_channel;i++)pa[0]+=pa[i];
	for(i=1;i<=n_channel;i++)pa[i]/=pa[0];
	for(i=1;i<=n_channel;i++)
		for(j=1;j<=i;j++)
			pi[i]+=pa[j];

	pi[0]=rand();
	pi[0]/=RAND_MAX;

	for(i=1;i<6;i++)
		if(pi[0]<=pi[i])break;

	p_channel[i]();
	stime+=1.0/pa[0];
}