source: trunk/MagicSoft/Simulation/Detector/Starfield/star.cxx@ 341

Last change on this file since 341 was 341, checked in by petry, 25 years ago
First version of the Starfield Generator in this repository. Fully functional version using the SKY2000 star catalog.
File size: 7.8 KB
Line 
1#include "star.hxx"
2
3star::star(){ // constructor (set invalid values)
4 icatnum = -999;
5 ra_h = -999.;
6 dec_deg = -999.;
7 umag = -999.;
8 bmag = -999.;
9 vmag = -999.;
10 rmag = -999.;
11 u = -999.;
12 v = -999.;
13 ra_rad = -999.;
14 dec_rad = -999.;
15}
16
17int star::readstar(FILE *fp, int verbose){ // read one line of the SKY2000 V2.0
18 // catalog and extract the interesting
19 // data
20 int ira_hours, ira_min;
21 int idec_degrees, idec_arcmin;
22 float ra_sec, dec_arcsec;
23 char catline[SKY2000LINELENGTH + 1];
24 char *pos;
25 char c2[3];
26 char c3[4];
27 char c6[7];
28 char c7[8];
29 char c8[9];
30
31 pos = catline;
32
33 if( fgets( pos , SKY2000LINELENGTH + 1, fp) == NULL ){
34 return(FALSE);
35 }
36
37 if(verbose > 2) fprintf(stdout, "%s\n", catline);
38
39 pos = catline + 27;
40 strncpy(c8, pos, 8);
41 sscanf(c8, "%d", &icatnum);
42
43 pos = catline + 118;
44 strncpy(c2, pos, 2);
45 sscanf(c2, "%d", &ira_hours);
46
47 pos = catline + 120;
48 strncpy(c2, pos, 2);
49 sscanf(c2, "%d", &ira_min);
50
51 pos = catline + 122;
52 strncpy(c7, pos, 7);
53 sscanf(c7, "%f", &ra_sec);
54
55 pos = catline + 129;
56 strncpy(c3, pos, 3);
57 if( c3[1] == ' ' ){
58 c3[1] = '0';
59 }
60 if( c3[2] == ' ' ){
61 c3[2] = '0';
62 }
63 sscanf(c3, "%d", &idec_degrees);
64
65 pos = catline + 132;
66 strncpy(c2, pos, 2);
67 sscanf(c2, "%d", &idec_arcmin);
68
69 pos = catline + 134;
70 strncpy(c6, pos, 6);
71 sscanf(c6, "%f", &dec_arcsec);
72
73 pos = catline + 231;
74 strncpy(c6, pos, 6);
75 sscanf(c6, "%f", &vmag);
76
77 pos = catline + 251;
78 strncpy(c6, pos, 6);
79 sscanf(c6, "%f", &bmag);
80
81 pos = catline + 271;
82 strncpy(c6, pos, 6);
83 sscanf(c6, "%f", &umag);
84
85 ra_h = ira_hours + ira_min/60. + ra_sec/3600.;
86 dec_deg = idec_degrees + idec_arcmin/60. + dec_arcsec/3600.;
87 ra_rad = ra_h * PI / 12.;
88 dec_rad = dec_deg * PI /180.;
89
90 if (verbose > 2) fprintf(stdout, "extracted: %d %d %d %f %d %d %f %f %f %f\n",
91 icatnum, ira_hours, ira_min, ra_sec,
92 idec_degrees, idec_arcmin, dec_arcsec,
93 umag, bmag, vmag);
94
95 return(TRUE);
96}
97
98int star::printstar(){ // write one star's parameters
99 fprintf(stdout, "%d %f %f %f %f %f %f\n",
100 icatnum, ra_h, dec_deg, umag, bmag, vmag, rmag);
101 return(0);
102}
103
104//----------------------------------------------------------------------------
105// @name calcmissingmags
106//
107// @desc calculate the magnitudes for those wavebands in which no data
108// @desc is available assuming a black body and using the V and B mags
109//
110//----------------------------------------------------------------------------
111
112float star::calcmissingmags(int verbose) { // returns effective temperature; -1. = not possible
113
114 float temp;
115 float tprime;
116 float nu1_Hz, nu2_Hz, bflux, vflux, rflux, xmag;
117
118
119
120 if(vmag > -100.){ // valid vmag available
121 if(bmag < -100.){ // no valid bmag available
122 cout << "Warning: star no. " << icatnum <<
123 " has no Bmag measurement. Using Bmag = Vmag = "<<vmag<<"\n";
124 bmag = vmag;
125 }
126 }
127 else{
128 cout << "Warning: star no. " << icatnum <<
129 " has no Vmag measurement.\n";
130 return(-1.);
131 }
132
133 // calculate the star temperature using approximation from
134 // Kitchin, C.R., Astrophysical Techniques, 2nd ed., equ. 3.1.24
135
136 if((bmag-vmag) > -0.2){
137 temp = 8540. / ( (bmag-vmag) + 0.865 );
138 if (verbose > 1) cout << "Star temperature from B-V: T = " << temp << "K\n";
139 }
140 else{
141 temp = 12000.;
142 if (verbose > 1) cout << "Star temperature from B-V: T > " << temp << "K\n";
143 }
144
145 // calculate an effective temperature for the Rmag calculation
146 // tprime = T * k / h
147
148 nu1_Hz = LIGHTSPEED_mps/((VLMIN_nm+VLMAX_nm)/2.*1e-9);
149 nu2_Hz = LIGHTSPEED_mps/((BLMAX_nm+BLMAX_nm)/2.*1e-9);
150 vflux = pow(10.,-0.4*vmag-22.42);
151 bflux = pow(10.,-0.4*bmag-22.42);
152 tprime = (nu2_Hz - nu1_Hz) / ( log(vflux/bflux) - 3. * log(nu1_Hz/nu2_Hz) );
153 if (verbose > 1) cout << "Tprime = " << tprime/1.38e-23*6.62e-34 << "\n";
154
155
156 if( umag < -100. ){ // umag could not be read
157 if (verbose) cout << "Warning: star no. " << icatnum <<
158 " has no Umag measurement. Calculating it from its Vmag = " << vmag << "\n";
159 if (verbose) cout << " and Bmag = " << bmag << " assuming standard colour-colour-plot ... ";
160
161 if((bmag-vmag) > 1.4){
162 umag = bmag * 0.9;
163 }
164 else{
165 if((bmag-vmag) > 0.5){
166 umag = -0.5 + 1.37 * (bmag - vmag) + bmag;
167 }
168 else{
169 if((bmag-vmag) <= 0.){
170 umag = 4.07 * (bmag - vmag) + bmag;
171 }
172 else{
173 umag = 0.175 * (bmag - vmag) + bmag;
174 }
175 }
176 }
177
178 if (verbose) cout << " result Umag = " << umag << "\n";
179
180 if( umag < 5.0 ){
181 cout << "Warning: star no. " << icatnum << " is bright (Vmag =" << vmag << ", Bmag = "
182 << bmag << ")\n and has no Umag measurement. Estimated Umag is "<< umag <<"\n";
183 }
184
185 }
186 else{ // umag available
187 if (verbose > 1) {
188 cout << "Test: star no. " << icatnum <<
189 " has Umag = " << umag <<". Calculating it from its Vmag = " << vmag << "\n";
190 cout << " and Bmag = " << bmag << " assuming standard colour-colour-plot ...\n ";
191
192 if((bmag-vmag) > 1.4){
193 xmag = bmag * 0.9;
194 }
195 else{
196 if((bmag-vmag) > 0.5){
197 xmag = -0.5 + 1.37 * (bmag - vmag) + bmag;
198 }
199 else{
200 if((bmag-vmag) <= 0.){
201 xmag = 4.07 * (bmag - vmag) + bmag;
202 }
203 else{
204 xmag = 0.175 * (bmag - vmag) + bmag;
205 }
206 }
207 }
208 if (verbose > 2)
209 cout << "TEST " << umag <<" "<< xmag << " " << temp << " " << bmag << " " << vmag << "\n";
210
211 cout << " result Umag = " << xmag << "\n";
212 }
213
214 }
215
216 if( rmag < -100. ){ // rmag not present (it's not part of the catalog)
217
218 rflux = vflux * pow((VLMIN_nm + VLMAX_nm)/(RLMIN_nm + RLMAX_nm),3.) *
219 (exp(nu1_Hz/tprime) - 1.) /
220 (exp(LIGHTSPEED_mps/((RLMIN_nm+RLMAX_nm)/2.*1e-9)/tprime) - 1.);
221
222 rmag = (log10(rflux) + 22.42)/(-0.4);
223
224 }
225
226 return(tprime);
227}
228
229
230//----------------------------------------------------------------------------
231// @name mag_nphot
232//
233// @desc translates magnitudes in number of photons, using log(flux)= -0.4*m+22.42
234//
235//----------------------------------------------------------------------------
236
237int star::mag_nphot(int np[4], float inttime_s, float radius_m, int verbose) {
238
239 float bflux, vflux, uflux, rflux;
240 float bintensity, vintensity, uintensity, rintensity;
241 float unu1_Hz, unu2_Hz, bnu1_Hz, bnu2_Hz, vnu1_Hz, vnu2_Hz, rnu1_Hz, rnu2_Hz;
242
243 // The flux is given in Watt/m2*Hz.
244
245 uflux = pow(10.,-0.4*umag-22.42);
246 bflux = pow(10.,-0.4*bmag-22.42);
247 vflux = pow(10.,-0.4*vmag-22.42);
248 rflux = pow(10.,-0.4*rmag-22.42);
249
250 if (verbose) cout<<"MAGS "<<umag<<" "<<bmag<<" "<<vmag<<" "<<rmag<<endl;
251
252 unu1_Hz = LIGHTSPEED_mps/(ULMIN_nm*1e-9);
253 unu2_Hz = LIGHTSPEED_mps/(ULMAX_nm*1e-9);
254 bnu1_Hz = LIGHTSPEED_mps/(BLMIN_nm*1e-9);
255 bnu2_Hz = LIGHTSPEED_mps/(BLMAX_nm*1e-9);
256 vnu1_Hz = LIGHTSPEED_mps/(VLMIN_nm*1e-9);
257 vnu2_Hz = LIGHTSPEED_mps/(VLMAX_nm*1e-9);
258 rnu1_Hz = LIGHTSPEED_mps/(RLMIN_nm*1e-9);
259 rnu2_Hz = LIGHTSPEED_mps/(RLMAX_nm*1e-9);
260
261 // The intensity is given in number_of_photons/sec*m2. We obtain this units
262 // because we multiply by the conversion factor 1Joule/s=h*c/((lambda1+lambda2)/2)
263
264
265 uintensity = uflux*(unu1_Hz-unu2_Hz) * (ULMIN_nm+ULMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
266 bintensity = bflux*(bnu1_Hz-bnu2_Hz) * (BLMIN_nm+BLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
267 vintensity = vflux*(vnu1_Hz-vnu2_Hz) * (VLMIN_nm+VLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
268 rintensity = rflux*(rnu1_Hz-rnu2_Hz) * (RLMIN_nm+RLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
269
270 np[0] = (int)(uintensity * radius_m * radius_m * PI * inttime_s);
271 np[1] = (int)(bintensity * radius_m * radius_m * PI * inttime_s);
272 np[2] = (int)(vintensity * radius_m * radius_m * PI * inttime_s);
273 np[3] = (int)(rintensity * radius_m * radius_m * PI * inttime_s);
274
275 if (verbose) cout<<"NPH "<<np[0]<<" "<<np[1]<<" "<<np[2]<<" " <<np[3]<<endl;
276
277 numphot = np[0] + np[1] + np[2] + np[3];
278
279 return(numphot);
280
281}
Note: See TracBrowser for help on using the repository browser.