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

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