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

Last change on this file since 348 was 342, checked in by petry, 25 years ago
Changed calculation of the R magnitude to use the temperature derived from B-V rather than that assuming a black body.
File size: 7.9 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 << "Blackbody T = " << 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 temp = temp * 1.38e-23 / 6.62e-34; // * k / h
219
220 rflux = vflux * pow((VLMIN_nm + VLMAX_nm)/(RLMIN_nm + RLMAX_nm),3.) *
221 (exp(nu1_Hz/temp) - 1.) /
222 (exp(LIGHTSPEED_mps/((RLMIN_nm+RLMAX_nm)/2.*1e-9)/temp) - 1.);
223
224 rmag = (log10(rflux) + 22.42)/(-0.4);
225
226 }
227
228 return(tprime);
229}
230
231
232//----------------------------------------------------------------------------
233// @name mag_nphot
234//
235// @desc translates magnitudes in number of photons, using log(flux)= -0.4*m+22.42
236//
237//----------------------------------------------------------------------------
238
239int star::mag_nphot(int np[4], float inttime_s, float radius_m, int verbose) {
240
241 float bflux, vflux, uflux, rflux;
242 float bintensity, vintensity, uintensity, rintensity;
243 float unu1_Hz, unu2_Hz, bnu1_Hz, bnu2_Hz, vnu1_Hz, vnu2_Hz, rnu1_Hz, rnu2_Hz;
244
245 // The flux is given in Watt/m2*Hz.
246
247 uflux = pow(10.,-0.4*umag-22.42);
248 bflux = pow(10.,-0.4*bmag-22.42);
249 vflux = pow(10.,-0.4*vmag-22.42);
250 rflux = pow(10.,-0.4*rmag-22.42);
251
252 if (verbose) cout<<"MAGS "<<umag<<" "<<bmag<<" "<<vmag<<" "<<rmag<<endl;
253
254 unu1_Hz = LIGHTSPEED_mps/(ULMIN_nm*1e-9);
255 unu2_Hz = LIGHTSPEED_mps/(ULMAX_nm*1e-9);
256 bnu1_Hz = LIGHTSPEED_mps/(BLMIN_nm*1e-9);
257 bnu2_Hz = LIGHTSPEED_mps/(BLMAX_nm*1e-9);
258 vnu1_Hz = LIGHTSPEED_mps/(VLMIN_nm*1e-9);
259 vnu2_Hz = LIGHTSPEED_mps/(VLMAX_nm*1e-9);
260 rnu1_Hz = LIGHTSPEED_mps/(RLMIN_nm*1e-9);
261 rnu2_Hz = LIGHTSPEED_mps/(RLMAX_nm*1e-9);
262
263 // The intensity is given in number_of_photons/sec*m2. We obtain this units
264 // because we multiply by the conversion factor 1Joule/s=h*c/((lambda1+lambda2)/2)
265
266
267 uintensity = uflux*(unu1_Hz-unu2_Hz) * (ULMIN_nm+ULMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
268 bintensity = bflux*(bnu1_Hz-bnu2_Hz) * (BLMIN_nm+BLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
269 vintensity = vflux*(vnu1_Hz-vnu2_Hz) * (VLMIN_nm+VLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
270 rintensity = rflux*(rnu1_Hz-rnu2_Hz) * (RLMIN_nm+RLMAX_nm)*1e-9/2. /(PLANCK_si*LIGHTSPEED_mps);
271
272 np[0] = (int)(uintensity * radius_m * radius_m * PI * inttime_s);
273 np[1] = (int)(bintensity * radius_m * radius_m * PI * inttime_s);
274 np[2] = (int)(vintensity * radius_m * radius_m * PI * inttime_s);
275 np[3] = (int)(rintensity * radius_m * radius_m * PI * inttime_s);
276
277 if (verbose) cout<<"NPH "<<np[0]<<" "<<np[1]<<" "<<np[2]<<" " <<np[3]<<endl;
278
279 numphot = np[0] + np[1] + np[2] + np[3];
280
281 return(numphot);
282
283}
Note: See TracBrowser for help on using the repository browser.