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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
/*
Copyright (C) 2004-2023 Sergey Koposov
Email: skoposov AT ed DOT ac DOT uk
This file is part of Q3C.
Q3C 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.
Q3C 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 Q3C; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#include <string.h>
static int64_t rand_state = 1;
static const int64_t c = 12345;
static const int64_t m = ((int64_t)1) << 31;
static const int64_t a = 1103515245;
int64_t get_rand()
{
rand_state = ( a * rand_state + c ) % m;
return rand_state;
}
/* Run as
* $ ./gen_data 1 100
* The first argument of the program is the random seed for the pseudorandom
* sequence.
* The second number is the number of objects outputted
* Additional flags
* --withpm output pms as well
* --pmscale= maximum allowed pm in mas/yr
* --epoch= epoch of coordinates
* --randomepoch assign random epoch
*
*/
/* The random number sequence was based on Knuth's theorem A (from
* his second book)
*/
int main(int argc, char *argv[])
{
const int nrabins = 36000;
const int ndecbins = 18000;
const int ntotbins = nrabins * ndecbins; /* 2^x*3^y*5^z */
double corrections[ndecbins], total = 0, pmra, pmdec, pmscale = 1;
int npoints;
bool random_epoch = false;
double epoch, cur_epoch;
int i, extraarg;
char parsing_error = 1;
bool withpm = false;
// first argument is the seed and then number of points to generate
if (argc >= 3)
{
rand_state = atoi(argv[1]);
get_rand(); // advance one step
npoints = atoi(argv[2]);
parsing_error = 0;
for (extraarg = 0; extraarg < (argc - 3); extraarg++)
{
char *curarg = argv[3 + extraarg];
if (strncmp(curarg,"--randomepoch", 13) == 0) {random_epoch = true;}
if (strncmp(curarg,"--withpm", 9) == 0) {withpm = true;}
if (strncmp(curarg,"--pmscale=", 10) == 0)
{
if (sscanf(curarg, "--pmscale=%lf", &pmscale) == 0)
{
fprintf(stderr, "Formatting error of pmscale\n");
exit(1);
}
}
if (strncmp(curarg,"--epoch=",8) == 0)
{
sscanf(curarg, "--epoch=%lf", &epoch);
}
}
}
if (parsing_error)
{
fprintf(stderr, "Wrong arguments!\n"
"MUST be ./gen_data [RANDOM SEED] [NPOINTS] [PROPERMOTIONSCALE(optional)]");
exit(1);
}
for (i = 0; i < ndecbins; i++)
/* weights in order to have cosine distribution of declinations
* corresponding to uniform distribution on the sky */
{
corrections[i] = cos((-90. + (180. / ndecbins) * (i + 0.5)) *
M_PI / 180.);
}
int npointsleft = npoints;
while (npointsleft)
{
int64_t ra = (int64_t)(get_rand() * 1. / m * nrabins);
int64_t dec = (int64_t)(get_rand() * 1. / m * ndecbins);
if (withpm )
{
pmra = ((get_rand() * 1. / m) * 2 - 1) * pmscale;
pmdec = ((get_rand() * 1. / m) * 2 - 1) * pmscale;
if (random_epoch)
{
cur_epoch = ((get_rand() * 1. / m) ) * 20 + 1980;
}
else
{
cur_epoch = epoch;
}
}
if (get_rand() < (corrections[dec] * m))
{
if (withpm)
{
printf("%f %f %f %f %f\n", ra * (360. / nrabins),
-90 + dec * (180. / ndecbins), pmra, pmdec, cur_epoch);
}
else
{
printf("%f %f\n", ra * (360. / nrabins),
-90 + dec * (180. / ndecbins));
}
npointsleft--;
}
}
}
|
