/
generate_catalog.c
98 lines (83 loc) · 3.1 KB
/
generate_catalog.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
/*
* Copyright (c) 2020 Marzia Rivi
*
* This file is part of RadioLensfit2.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_rng.h>
#include "datatype.h"
#include "utils.h"
#include "default_params.h"
#include "distributions.h"
#include "generate_random_values.h"
#include "generate_catalog.h"
// Generate galaxy catalog ordered by source flux
// 2NP = number of sampled orientations (points on the circle of radius |e|) for each ellipticity module
unsigned int galaxy_catalog(unsigned int nge, int NP, double fov_eff, double Rmin, double Rmax, double Fmin, double Fmax, double* gflux, double* gscale, double* ge1, double *ge2, double *l, double *m)
{
//setup random number generator
const gsl_rng_type * G;
gsl_rng * gen;
G = gsl_rng_mt19937; // Mersenne Twister
gen = gsl_rng_alloc(G);
unsigned long int seed = random_seed();
gsl_rng_set(gen,seed);
unsigned int my_gal = nge;
unsigned int diffgal = my_gal/(2*NP);
my_gal = diffgal*2*NP;
// generate flux values
double *gflux2 = new double[diffgal];
generate_random_data(gen,diffgal,gflux2,Fmin,Fmax,flux_CDF,M_EXP);
// sort flux values, so that to generate a population ordered by flux and therefore fitting sources by decreasing flux order
gsl_sort(gflux2,1,diffgal); // sorting ascending order
unsigned int gal = 0;
for (unsigned int i=0; i<diffgal; i++)
for (int k = 0; k < 2*NP; k++)
{
gflux[gal] = gflux2[diffgal-i-1];
gal++;
}
delete[] gflux2;
// generate scalelength
double mu, scalelength;
for (unsigned int g = 0; g < my_gal; g++)
{
if (g%(2*NP) == 0)
{
mu = scale_mean(gflux[g]); //power law relation between flux and scalelength
generate_random_data(gen,1,&scalelength,Rmin,Rmax,r_CDF,mu);
}
gscale[g] = scalelength;
}
// generate ellipticities
generate_ellipticity(gen,diffgal,NP,ge1,ge2);
// generate positions
// uniformly random positions in RAD in a disk area
// http://mathworld.wolfram.com/DiskPointPicking.html
double radius,orient;
for (unsigned int gal=0; gal<my_gal; gal++)
{
radius = sqrt(gsl_rng_uniform(gen))*0.5*fov_eff;
orient = gsl_rng_uniform(gen)*2*PI;
l[gal] = radius*cos(orient);
m[gal] = radius*sin(orient);
}
gsl_rng_free(gen);
return my_gal;
}