создание Powerlaw с использованием Монте-Карло в C ++ - PullRequest
Новая
       95

создание Powerlaw с использованием Монте-Карло в C ++

0 голосов
/ 08 ноября 2019

В настоящее время я пытаюсь создать распределение powerlaw + gaussian с использованием кода Монте-Карло, который я написал: 

#include <cstdlib>
#include <stdio.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>   // format manipulation
#include <string>
#include <vector>
#include <algorithm>
#include <sstream>
#include <ctime>


#define PI 3.14159265

using namespace std;

int main(int argc, char *argv[]) {

double ts=108*24*3600;//tstart obs 
double T=50000;//obs duration

double r; 
double A=3.15;
double B=-0.47; 
double rstar=ts+T*2/3; 
double Rmax; 
int counts=100000; 
double C=0.005; 
double sigma= 1500;

double maxpow=A*pow(ts,B)+C/((sigma)*sqrt(2*PI))*exp(-2*pow(((ts-rstar)/sigma),2));//ok 
double maxflare=A*pow(rstar,B)+C/((sigma)*sqrt(2*PI));//ok 
double minpow=A*pow(ts+T,B)+C/((sigma)*sqrt(2*PI))*exp(-2*pow(((ts+T-rstar)/sigma),2));//ok 
double powgaussian;

for(int i=0;i<counts;i++){

extraction: //random time in the observation 
r = ((double) rand() / (RAND_MAX))*T+ts;//ok

//forma LC 
powgaussian = A*pow(r,B)+ C/((sigma)*sqrt(2*PI))  *  exp(-2*pow(  ((r-rstar)/sigma),2) );//ok

//max della LC  
if ( maxpow >= maxflare )   Rmax=((double) rand() / (RAND_MAX))*(maxpow-minpow)+minpow;//ok
else if (maxpow < maxflare) Rmax=((double) rand() / (RAND_MAX))*(maxflare-minpow)+minpow;//ok


if (Rmax>powgaussian)  goto extraction;

ofstream out("test.dat", ios::app); out << r << "   "<<powgaussian<<endl; out.close(); //ok
}

в выходном файле, если я попытаюсь построить извлеченный «r» против «powgaussian»(r) "значение, которое я получаю именно так, как я ожидаю the powgaussian variable is calculated correctly

, однако, если я создаю гистограмму извлеченных значений" r ", они не следуют" паугауссовой "формевообще the distribution of extracted

Я не понимаю причину этого несоответствия, мой код кажется правильным, и у меня заканчиваются варианты. Кто-нибудь имеет представление о том, что я делаю неправильно? (log-log оси на обоих изображениях)

Спасибо

...