Line data Source code
1 : //# FluxCalc_SS_JPL_Butler.cc: Implementation of FluxCalc_SS_JPL_Butler.h
2 : //# Copyright (C) 2010
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //----------------------------------------------------------------------------
27 :
28 : //#include <components/ComponentModels/FluxStandard.h>
29 : #include <components/ComponentModels/FluxCalc_SS_JPL_Butler.h>
30 : #include <components/ComponentModels/ComponentType.h>
31 : #include <casacore/casa/Containers/Record.h>
32 : #include <casacore/casa/BasicMath/Math.h>
33 : #include <casacore/casa/BasicSL/String.h>
34 : #include <casacore/casa/Quanta.h>
35 : #include <casacore/casa/Logging/LogIO.h>
36 : #include <casacore/casa/OS/Directory.h>
37 : #include <casacore/casa/OS/DirectoryIterator.h>
38 : #include <casacore/casa/System/Aipsrc.h>
39 : #include <sstream>
40 : #include <iomanip>
41 : #include <casacore/measures/Measures.h>
42 : #include <casacore/measures/Measures/MEpoch.h>
43 : #include <casacore/measures/Measures/MCEpoch.h>
44 : #include <casacore/measures/Measures/MDirection.h>
45 : #include <casacore/measures/Measures/MFrequency.h>
46 : #include <casacore/measures/Measures/MeasComet.h>
47 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
48 : #include <casacore/tables/Tables/Table.h>
49 : #include <casacore/tables/Tables/TableRecord.h>
50 : #include <casacore/tables/Tables/ScalarColumn.h>
51 : #include <casatools/Config/State.h>
52 :
53 : using namespace casacore;
54 : namespace casa { //# NAMESPACE CASA - BEGIN
55 :
56 : // Recommended constructor.
57 7 : FluxCalc_SS_JPL_Butler::FluxCalc_SS_JPL_Butler(const String& objname,
58 7 : const MEpoch& time) :
59 : name_p(objname),
60 : hasName_p(true),
61 : time_p(time),
62 : hasTime_p(true),
63 : hasEphemInfo_p(false),
64 : hertz_p("Hz"),
65 : has_ra_p(false),
66 : has_dec_p(false),
67 7 : has_illu_p(false)
68 : {
69 7 : hasObjNum_p = setObjNum();
70 7 : }
71 :
72 0 : FluxCalc_SS_JPL_Butler::FluxCalc_SS_JPL_Butler() :
73 : name_p(""),
74 : hasName_p(false),
75 : time_p(MEpoch()),
76 : hasTime_p(false),
77 : hasEphemInfo_p(false),
78 0 : hertz_p("Hz")
79 : {
80 : // Default constructor for making arrays, etc.
81 0 : hasObjNum_p = false;
82 0 : objnum_p = FluxCalc_SS_JPL_Butler::N_KNOWN;
83 0 : }
84 :
85 7 : FluxCalc_SS_JPL_Butler::~FluxCalc_SS_JPL_Butler()
86 : {
87 : // Default destructor
88 7 : }
89 :
90 7 : Bool FluxCalc_SS_JPL_Butler::setObjNum()
91 : {
92 21 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "setObjNum"));
93 :
94 7 : if(!hasName_p){
95 : os << LogIO::SEVERE
96 : << "Please provide the source name."
97 0 : << LogIO::POST;
98 0 : return false;
99 : }
100 :
101 7 : String lname = name_p;
102 7 : lname.downcase();
103 7 : Bool matched = true;
104 :
105 7 : if(lname == "mercury")
106 0 : objnum_p = FluxCalc_SS_JPL_Butler::Mercury;
107 7 : else if(lname == "venus")
108 0 : objnum_p = FluxCalc_SS_JPL_Butler::Venus;
109 7 : else if(lname == "mars")
110 1 : objnum_p = FluxCalc_SS_JPL_Butler::Mars;
111 6 : else if(lname == "jupiter")
112 0 : objnum_p = FluxCalc_SS_JPL_Butler::Jupiter;
113 6 : else if(lname == "io")
114 0 : objnum_p = FluxCalc_SS_JPL_Butler::Io;
115 6 : else if(lname == "ganymede")
116 0 : objnum_p = FluxCalc_SS_JPL_Butler::Ganymede;
117 6 : else if(lname == "europa")
118 0 : objnum_p = FluxCalc_SS_JPL_Butler::Europa;
119 6 : else if(lname == "callisto")
120 0 : objnum_p = FluxCalc_SS_JPL_Butler::Callisto;
121 6 : else if(lname == "uranus")
122 0 : objnum_p = FluxCalc_SS_JPL_Butler::Uranus;
123 6 : else if(lname == "neptune")
124 0 : objnum_p = FluxCalc_SS_JPL_Butler::Neptune;
125 6 : else if(lname == "triton")
126 0 : objnum_p = FluxCalc_SS_JPL_Butler::Triton;
127 6 : else if(lname == "pluto")
128 0 : objnum_p = FluxCalc_SS_JPL_Butler::Pluto;
129 6 : else if(lname == "titan")
130 5 : objnum_p = FluxCalc_SS_JPL_Butler::Titan;
131 1 : else if(lname == "ceres")
132 0 : objnum_p = FluxCalc_SS_JPL_Butler::Ceres;
133 1 : else if(lname == "pallas")
134 0 : objnum_p = FluxCalc_SS_JPL_Butler::Pallas;
135 1 : else if(lname == "vesta")
136 0 : objnum_p = FluxCalc_SS_JPL_Butler::Vesta;
137 1 : else if(lname == "juno")
138 0 : objnum_p = FluxCalc_SS_JPL_Butler::Juno;
139 1 : else if(lname == "victoria")
140 0 : objnum_p = FluxCalc_SS_JPL_Butler::Victoria;
141 1 : else if(lname == "davida")
142 0 : objnum_p = FluxCalc_SS_JPL_Butler::Davida;
143 : else{
144 : os << LogIO::SEVERE
145 1 : << "Sorry, no flux density model for " << name_p
146 : << "\n (not even a rudimentary one)."
147 1 : << LogIO::POST;
148 1 : matched = false;
149 : }
150 :
151 : // Changing the object invalidates the cached ephemeris info (if any).
152 7 : hasEphemInfo_p = false;
153 :
154 7 : return matched;
155 : }
156 :
157 0 : Bool FluxCalc_SS_JPL_Butler::getName(String& output) const
158 : {
159 0 : if(!hasName_p)
160 0 : return false;
161 0 : output = name_p;
162 0 : return true;
163 : }
164 :
165 0 : Bool FluxCalc_SS_JPL_Butler::getTime(MEpoch& output) const
166 : {
167 0 : if(!hasTime_p)
168 0 : return false;
169 0 : output = time_p;
170 0 : return true;
171 : }
172 :
173 0 : void FluxCalc_SS_JPL_Butler::setTime(const MEpoch& time)
174 : {
175 0 : time_p = time;
176 :
177 : // Changing the time *possibly* invalidates the cached ephemeris info (if
178 : // any). Leave it up to readEphem() to decide whether or not the ephemeris
179 : // info is valid (i.e. the new time is close enough to the old time).
180 0 : hasEphemInfo_p = false;
181 0 : }
182 :
183 17 : ComponentType::Shape FluxCalc_SS_JPL_Butler::getShape(Double& angdiam)
184 : {
185 17 : if(!hasEphemInfo_p && !readEphem())
186 3 : return ComponentType::UNKNOWN_SHAPE;
187 :
188 14 : angdiam = 2.0 * mean_rad_p / delta_p;
189 14 : return ComponentType::DISK;
190 : }
191 :
192 0 : MDirection FluxCalc_SS_JPL_Butler::getDirection()
193 : {
194 0 : if((!hasEphemInfo_p && !readEphem()) || !(has_ra_p && has_dec_p))
195 0 : return MDirection();
196 :
197 0 : return MDirection(MVDirection(Quantity(ra_p, "deg"),
198 0 : Quantity(dec_p, "deg")), MDirection::J2000);
199 : }
200 :
201 0 : Double FluxCalc_SS_JPL_Butler::getHeliocentricDist()
202 : {
203 : Double dist;
204 :
205 0 : if(!hasEphemInfo_p && !readEphem())
206 0 : dist = -1.0;
207 0 : if(!has_r_p)
208 0 : dist = -1.0;
209 : else
210 0 : dist = r_p;
211 0 : return dist;
212 : }
213 :
214 0 : uInt FluxCalc_SS_JPL_Butler::n_known() const
215 : {
216 0 : return N_KNOWN;
217 : }
218 :
219 7 : Bool FluxCalc_SS_JPL_Butler::readEphem()
220 : {
221 21 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "readEphem"));
222 :
223 7 : if(!hasName_p || !hasTime_p){
224 : os << LogIO::SEVERE
225 : << "The source and time have not been set."
226 0 : << LogIO::POST;
227 0 : return false;
228 : }
229 :
230 : // Try to find a matching JPL-Horizons ephemeris table.
231 : // Note: these are not the same as the DE200 and DE405 JPL tables used
232 : // by measures to get the direction to a planet.
233 : // There may be more than one because of overlapping date ranges.
234 7 : const String tabpat(Regex::makeCaseInsensitive(name_p +
235 14 : "_[-0-9.]+-[-0-9.]+[ydhms].+\\.tab"));
236 :
237 : // Look for tabpat, in order, in '.', user.ephemerides.directory, and
238 : // Aipsrc::findDir(horpath, "data/ephemerides/JPL-Horizons").
239 7 : uInt nephdirs = 1;
240 14 : String userpath;
241 7 : Bool foundUserDir = Aipsrc::find(userpath, "user.ephemerides.directory");
242 7 : if(foundUserDir)
243 0 : ++nephdirs;
244 14 : String resolvepath;
245 7 : Bool foundResolveDir = false;
246 7 : resolvepath = casatools::get_state( ).resolve("ephemerides/JPL-Horizons");
247 7 : if (resolvepath != "ephemerides/JPL-Horizons") {
248 7 : foundResolveDir = true;
249 7 : ++nephdirs;
250 : }
251 14 : String horpath;
252 7 : Bool foundStd = Aipsrc::findDir(horpath, "data/ephemerides/JPL-Horizons");
253 7 : if(foundStd)
254 0 : ++nephdirs;
255 :
256 7 : int nephindex = 0;
257 14 : Vector<String> ephdirs(nephdirs);
258 7 : ephdirs[nephindex++] = ".";
259 7 : if(foundUserDir)
260 0 : ephdirs[nephindex++] = userpath;
261 7 : if(foundResolveDir)
262 7 : ephdirs[nephindex++] = resolvepath;
263 7 : if(foundStd)
264 0 : ephdirs[nephindex++] = horpath;
265 :
266 7 : Bool foundObj = false;
267 7 : Bool found = false; // = foundObj && right time.
268 14 : Path path;
269 14 : String edir(".");
270 21 : for(uInt pathnum = 0; pathnum < nephdirs && !found; ++pathnum){
271 14 : edir = ephdirs[pathnum];
272 :
273 : os << LogIO::NORMAL2
274 : << "Looking for an ephemeris table matching " << tabpat
275 : << "\n\tin " << edir
276 14 : << LogIO::POST;
277 :
278 28 : Directory hordir(edir);
279 28 : DirectoryIterator dirIter(hordir, Regex(tabpat));
280 14 : uInt firstTimeStart = name_p.length() + 1; // The + 1 is for the _.
281 42 : Regex timeUnitPat("[ydhms]");
282 :
283 35 : while(!dirIter.pastEnd()){
284 25 : path = dirIter.name();
285 25 : foundObj = true;
286 25 : String basename(path.baseName());
287 :
288 : // Look for, respectively, the positions of '--', 'd', and '.' in
289 : // '-12345--67890dUTC.tab'. Note that, just to be tricky, the times in
290 : // this example are negative.
291 25 : uInt firstTimeLen = basename.find('-', firstTimeStart + 1) - firstTimeStart;
292 50 : uInt lastTimeLen = basename.find(timeUnitPat,
293 25 : firstTimeStart + firstTimeLen + 1)
294 25 : - firstTimeStart - firstTimeLen - 1;
295 25 : uInt unitPos = firstTimeStart + firstTimeLen + 1 + lastTimeLen;
296 :
297 25 : Double firstTime = String::toDouble(basename.substr(firstTimeStart, firstTimeLen));
298 25 : Double lastTime = String::toDouble(basename.substr(firstTimeStart + firstTimeLen + 1,
299 : lastTimeLen));
300 25 : Unit unit(basename[unitPos]);
301 25 : String ref(basename.substr(unitPos + 1,
302 25 : basename.find('.', unitPos + 1) - unitPos - 1));
303 :
304 : os << LogIO::DEBUG1
305 : << basename << ": (first, last)time = ("
306 : << firstTime << ", " << lastTime << ")" << unit.getName()
307 : << " " << ref
308 25 : << LogIO::POST;
309 :
310 : MEpoch::Types refEnum;
311 25 : Bool refIsValid = MEpoch::getType(refEnum, ref);
312 :
313 25 : if(refIsValid){
314 15 : MEpoch::Convert mtimeToDirFrame(time_p, MEpoch::Ref(refEnum));
315 15 : MEpoch mtimeInDirFrame(mtimeToDirFrame());
316 15 : Double dtime = mtimeInDirFrame.get(unit).getValue();
317 :
318 15 : if(dtime <= lastTime && dtime >= firstTime){
319 4 : found = true;
320 4 : break;
321 : }
322 : }
323 : // else maybe tabpat isn't specific enough. Don't panic yet.
324 :
325 21 : ++dirIter;
326 : }
327 : }
328 :
329 7 : if(!found){
330 3 : os << LogIO::SEVERE;
331 3 : if(foundObj)
332 3 : os << "Found an ephemeris for " << name_p << ", but not";
333 : else
334 0 : os << "Could not find an ephemeris table for " << name_p;
335 :
336 : // MEpoch cannot directly << to a LogIO.
337 3 : os << " at ";
338 3 : os.output() << MEpoch::Convert(time_p, MEpoch::Ref(MEpoch::UTC))();
339 3 : os << LogIO::POST;
340 3 : return false;
341 : }
342 : else{
343 : os << LogIO::NORMAL
344 4 : << "Using ephemeris table " << path.baseName()
345 4 : << LogIO::POST;
346 : }
347 :
348 : // path.absoluteName() is liable to give something like cwd +
349 : // path.baseName(), because path was never given horpath.
350 12 : const String abspath(edir + "/" + path.baseName());
351 :
352 4 : if(!Table::isReadable(abspath)){
353 : os << LogIO::SEVERE
354 : << abspath << " is not a readable table."
355 0 : << LogIO::POST;
356 0 : return false;
357 : }
358 :
359 8 : const Table tab(abspath);
360 8 : const TableRecord ks(tab.keywordSet());
361 :
362 4 : Bool got_q = true;
363 4 : temperature_p = MeasComet::get_Quantity_keyword(ks, "T_mean", "K", got_q);
364 4 : if(!got_q)
365 4 : temperature_p = -1; // Hopefully a model for the obj will supply a
366 : // temperature later.
367 4 : mean_rad_p = MeasComet::get_Quantity_keyword(ks, "meanrad", "AU", got_q);
368 4 : if(!got_q){
369 0 : mean_rad_p = -1.0;
370 : os << LogIO::SEVERE // Remove/modify this when it starts supporting triaxiality.
371 : << "The table is missing the meanrad keyword, needed to calculate the apparent diameter."
372 0 : << LogIO::POST;
373 0 : return false;
374 : }
375 :
376 : // Find the row numbers with the right MJDs.
377 8 : ScalarColumn<Double> mjd(tab, "MJD");
378 : uInt rowbef;
379 : uInt rowclosest;
380 : uInt rowaft;
381 4 : if(!get_row_numbers(rowbef, rowclosest, rowaft, mjd)){
382 : os << LogIO::SEVERE
383 : << "The table does not appear to cover the right time."
384 0 : << LogIO::POST;
385 0 : return false;
386 : }
387 :
388 4 : Double tm1 = mjd(rowbef);
389 4 : Double t0 = mjd(rowclosest);
390 4 : Double tp1 = mjd(rowaft);
391 4 : Double f = time_p.get("d").getValue() - t0;
392 4 : Double dt = tp1 - tm1;
393 4 : Double tp1mt0 = tp1 - t0;
394 4 : Double t0mtm1 = t0 - tm1;
395 :
396 : // The distance from Earth to the object, in AU, is mandatory.
397 : // JPL calls it delta, and MeasComet calls it Rho.
398 4 : hasEphemInfo_p = found && get_interpolated_value(delta_p, "Rho",
399 : tab, rowbef, rowclosest,
400 : rowaft, f, dt, tp1mt0,
401 : t0mtm1, true);
402 :
403 : // Heliocentric distance, in AU.
404 4 : has_r_p = get_interpolated_value(r_p, "r", tab, rowbef, rowclosest, rowaft,
405 : f, dt, tp1mt0, t0mtm1, false);
406 :
407 : // Illumination, in %.
408 4 : has_illu_p = get_interpolated_value(illu_p, "illu", tab, rowbef, rowclosest,
409 : rowaft, f, dt, tp1mt0, t0mtm1, false);
410 4 : if(has_illu_p)
411 4 : has_illu_p *= 0.01; // Convert it to a fraction.
412 :
413 : // RA, in deg.
414 4 : has_ra_p = get_interpolated_value(ra_p, "RA", tab, rowbef, rowclosest,
415 : rowaft, f, dt, tp1mt0, t0mtm1, false);
416 :
417 : // Declination, in deg.
418 4 : has_dec_p = get_interpolated_value(dec_p, "DEC", tab, rowbef, rowclosest,
419 : rowaft, f, dt, tp1mt0, t0mtm1, false);
420 :
421 4 : return found;
422 : }
423 :
424 20 : Bool FluxCalc_SS_JPL_Butler::get_interpolated_value(Double& val,
425 : const String& colname,
426 : const Table& tab,
427 : const uInt rowbef,
428 : const uInt rowclosest,
429 : const uInt rowaft,
430 : const Double f,
431 : const Double dt,
432 : const Double tp1mt0,
433 : const Double t0mtm1,
434 : const Bool verbose)
435 : {
436 20 : Bool foundIt = false;
437 40 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "get_interpolated_value"));
438 :
439 20 : if(tab.actualTableDesc().isColumn(colname)){
440 20 : Double myf = f;
441 20 : Double d2y = 0.0;
442 :
443 20 : ScalarColumn<Double> col(tab, colname);
444 20 : Double col_m1 = col(rowbef);
445 20 : Double col_0 = col(rowclosest);
446 20 : Double col_p1 = col(rowaft);
447 :
448 20 : if(dt > 0){
449 20 : myf /= dt;
450 20 : if(t0mtm1 > 0.0 && tp1mt0 > 0.0){
451 20 : d2y = (col_p1 - col_0) / tp1mt0;
452 20 : d2y -= (col_0 - col_m1) / t0mtm1;
453 20 : d2y *= dt;
454 : }
455 : }
456 : else{
457 0 : if(verbose){
458 : os << LogIO::NORMAL
459 : << "The table is not long enough for quadratic interpolation.\n"
460 : << "Nearest neighbor will be used."
461 0 : << LogIO::POST;
462 : }
463 0 : myf = 0.0;
464 : }
465 20 : val = col_0 + myf * (col_p1 - col_m1 + myf * d2y);
466 20 : foundIt = true;
467 : }
468 : else
469 : os << LogIO::NORMAL
470 : << "The table does not have a " << colname << " column."
471 0 : << LogIO::POST;
472 40 : return foundIt;
473 : }
474 :
475 4 : Bool FluxCalc_SS_JPL_Butler::get_row_numbers(uInt& rowbef, uInt& rowclosest,
476 : uInt& rowaft,
477 : const ScalarColumn<Double>& mjd)
478 : {
479 : // MeasComet requires a constant time increment, but since
480 : // FluxCalc_SS_JPL_Butler is expected to only need to use the time once, it's
481 : // not too expensive to allow tables with varying time increments. As long
482 : // as mjd is monotonically increasing, the search is at worst O(log(n)).
483 :
484 4 : Double mjd0 = mjd(0);
485 4 : Double dmjd = mjd0;
486 4 : Int ndates = mjd.nrow(); // Don't bother trying uInts in this
487 4 : Int step = 1; // function - it just leads to several
488 4 : Long rn = 0; // compiler warnings.
489 :
490 4 : Int ub = ndates - 1;
491 4 : Double the_time = time_p.get("d").getValue();
492 :
493 4 : if(mjd(ub) < the_time){
494 0 : return false;
495 : }
496 4 : else if(mjd(ub) == the_time){
497 0 : rn = ub;
498 0 : step = 0; // Prevents going through the while loop below.
499 : }
500 4 : Int lb = 0;
501 4 : if(mjd(0) > the_time){
502 0 : return false;
503 : }
504 4 : else if(mjd(0) == the_time){
505 0 : rn = 0;
506 0 : step = 0; // Prevents going through the while loop below.
507 : }
508 :
509 : Int i;
510 8 : for(i = 1; dmjd == mjd0 && i < ndates; ++i)
511 4 : dmjd = mjd(i);
512 4 : if(i > 1)
513 4 : --i;
514 4 : dmjd = (dmjd - mjd0) / i;
515 :
516 4 : if(dmjd > 0.0 && step){
517 4 : rn = lrint((the_time - mjd0) / dmjd);
518 4 : if(rn < 0)
519 0 : rn = 0;
520 4 : else if(rn > ndates)
521 0 : rn = ndates - 1;
522 : }
523 :
524 4 : Double mjdrn = mjd(rn);
525 4 : Bool increasing = mjdrn < the_time;
526 4 : Int paranoia = 0;
527 :
528 16 : while(step && paranoia < ndates){
529 12 : if(mjdrn < the_time){
530 4 : if(rn > lb)
531 4 : lb = rn;
532 4 : if(increasing){
533 4 : step *= 2;
534 : }
535 : else{
536 0 : step /= 2;
537 0 : increasing = true;
538 : }
539 : }
540 : else{
541 8 : if(rn < ub)
542 8 : ub = rn;
543 8 : if(increasing){
544 4 : step /= 2;
545 4 : increasing = false;
546 : }
547 : else{
548 4 : step *= 2;
549 : }
550 : }
551 12 : if(increasing){
552 4 : if(rn + step > ub)
553 0 : step = ub - rn - 1;
554 4 : rn += step;
555 : }
556 : else{
557 8 : if(rn - step < lb)
558 4 : step = rn - lb - 1;
559 8 : rn -= step;
560 : }
561 12 : mjdrn = mjd(rn);
562 12 : ++paranoia;
563 : }
564 4 : if(paranoia == ndates)
565 0 : return false;
566 :
567 4 : rowclosest = rn;
568 4 : rowaft = (rn < ndates - 1) ? rn + 1 : rn;
569 4 : rowbef = (rn > 0) ? rn - 1 : rn;
570 4 : return true;
571 : }
572 :
573 17 : ComponentType::Shape FluxCalc_SS_JPL_Butler::compute(Vector<Flux<Double> >& values,
574 : Vector<Flux<Double> >& errors,
575 : Double& angdiam,
576 : const Vector<MFrequency>& mfreqs,
577 : const Bool report)
578 : {
579 : // LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute"));
580 :
581 : // Calls readEphem() if necessary.
582 17 : ComponentType::Shape rettype(getShape(angdiam));
583 17 : if(rettype == ComponentType::UNKNOWN_SHAPE)
584 3 : return rettype;
585 :
586 14 : if(!hasObjNum_p){
587 0 : hasObjNum_p = setObjNum(); // Also has its own errmsgs.
588 0 : if(!hasObjNum_p)
589 0 : return ComponentType::UNKNOWN_SHAPE;
590 : }
591 :
592 14 : switch(objnum_p){
593 0 : case FluxCalc_SS_JPL_Butler::Venus:
594 0 : compute_venus(values, errors, angdiam, mfreqs);
595 0 : break;
596 0 : case FluxCalc_SS_JPL_Butler::Jupiter:
597 0 : compute_jupiter(values, errors, angdiam, mfreqs);
598 0 : break;
599 0 : case FluxCalc_SS_JPL_Butler::Uranus:
600 0 : compute_uranus(values, errors, angdiam, mfreqs);
601 0 : break;
602 0 : case FluxCalc_SS_JPL_Butler::Neptune:
603 0 : compute_neptune(values, errors, angdiam, mfreqs);
604 0 : break;
605 14 : default:
606 14 : Bool success = compute_constant_temperature(values, errors, angdiam, mfreqs,
607 : report);
608 14 : if(!success)
609 0 : return ComponentType::UNKNOWN_SHAPE;
610 : };
611 :
612 14 : return rettype;
613 : }
614 :
615 14 : void FluxCalc_SS_JPL_Butler::compute_BB(Vector<Flux<Double> >& values,
616 : Vector<Flux<Double> >& errors,
617 : const Double angdiam,
618 : const Vector<MFrequency>& mfreqs)
619 : {
620 14 : const uInt nfreqs = mfreqs.nelements();
621 28 : Quantum<Double> temperature(temperature_p, "K");
622 :
623 : // The real peak frequency is about 2.82 x this.
624 28 : Quantum<Double> freq_peak(QC::k( ) * temperature / QC::h( ));
625 :
626 28 : Quantum<Double> rocd2(0.5 * angdiam); // Dimensionless for now.
627 :
628 14 : rocd2 /= QC::c( ); // Don't put this in the c'tor, it'll give the wrong answer.
629 14 : rocd2 *= rocd2;
630 :
631 : // Frequency independent factor.
632 42 : Quantum<Double> freq_ind_fac(2.0e26 * QC::h( ) * C::pi * rocd2);
633 :
634 42 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_BB"));
635 : os << LogIO::DEBUG1
636 : << "angdiam = " << angdiam << " rad"
637 28 : << "\nrocd2 = " << rocd2.getValue() << rocd2.getUnit()
638 28 : << "\nfreq_ind_fac = " << freq_ind_fac.getValue() << freq_ind_fac.getUnit()
639 28 : << "\npeak freq = " << 2.82e-12 * freq_peak.get(hertz_p).getValue() << " THz"
640 : << "\ntemperature_p = " << temperature_p << " K"
641 14 : << "\nvalues[0].unit() = " << values[0].unit().getName()
642 : << "\nhertz_p = " << hertz_p.getName()
643 42 : << LogIO::POST;
644 :
645 28 : const Unit jy("Jy");
646 :
647 30740 : for(uInt f = 0; f < nfreqs; ++f){
648 61452 : Quantum<Double> freq(mfreqs[f].get(hertz_p));
649 :
650 30726 : values[f].setUnit(jy);
651 61452 : Double fd = (freq_ind_fac * freq * freq * freq).getValue() /
652 30726 : (exp((freq / freq_peak).getValue()) - 1.0);
653 : os << LogIO::DEBUG1
654 61452 : << "f.d.(" << freq.getValue() << " Hz" << ") = " << fd
655 30726 : << LogIO::POST;
656 30726 : values[f].setValue(fd);
657 30726 : errors[f].setValue(0.0);
658 : }
659 14 : }
660 :
661 0 : void FluxCalc_SS_JPL_Butler::compute_GB(Vector<Flux<Double> >& values,
662 : Vector<Flux<Double> >& errors,
663 : const Double angdiam,
664 : const Vector<MFrequency>& mfreqs,
665 : const Vector<Double>& temps)
666 : {
667 0 : const uInt nfreqs = mfreqs.nelements();
668 0 : Quantum<Double> rocd2(0.5 * angdiam); // Dimensionless for now.
669 :
670 0 : rocd2 /= QC::c( ); // Don't put this in the c'tor, it'll give the wrong answer.
671 0 : rocd2 *= rocd2;
672 :
673 : // Frequency independent factor.
674 0 : Quantum<Double> freq_ind_fac(2.0e26 * QC::h( ) * C::pi * rocd2);
675 :
676 0 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_GB"));
677 : os << LogIO::DEBUG1
678 : << "angdiam = " << angdiam << " rad"
679 0 : << "\nrocd2 = " << rocd2.getValue() << rocd2.getUnit()
680 0 : << "\nfreq_ind_fac = " << freq_ind_fac.getValue() << freq_ind_fac.getUnit()
681 0 : << LogIO::POST;
682 :
683 0 : const Unit jy("Jy");
684 :
685 0 : for(uInt f = 0; f < nfreqs; ++f){
686 0 : Quantum<Double> freq(mfreqs[f].get(hertz_p));
687 :
688 : // Guard against wayward extrapolations (possible with compute_venus()),
689 : // but do not emit a warning here; there may be many frequencies with bad
690 : // temperatures, and the warning should come from the calling function
691 : // anyway since it knows the limits of the model.
692 : //
693 : // I am not really claiming that the CMB temperature is the most reasonable
694 : // minimum temperature, but it is *a* reasonable mininum temperature, and I
695 : // want to avoid division by zero.
696 : //
697 0 : Quantum<Double> temperature(max(temps[f], 2.7), "K");
698 :
699 : // The real peak frequency is about 2.82 x this.
700 0 : Quantum<Double> freq_peak(QC::k( ) * temperature / QC::h( ));
701 :
702 0 : values[f].setUnit(jy);
703 0 : Double fd = (freq_ind_fac * freq * freq * freq).getValue() /
704 0 : (exp((freq / freq_peak).getValue()) - 1.0);
705 0 : values[f].setValue(fd);
706 0 : errors[f].setValue(0.0);
707 :
708 : // Take this out when it's served its purpose, since it's in a possibly
709 : // long loop.
710 : // os << LogIO::DEBUG2
711 : // << "f = 0 (" << 1e-12 * freq.get(hertz_p).getValue() << " THz):\n"
712 : // << "temperature = " << temps[f] << " K\n"
713 : // << "freq_peak = " << 1e-12 * freq_peak.get(hertz_p).getValue() << " THz\n"
714 : // << "f.d. = " << fd
715 : // << LogIO::POST;
716 : }
717 : os << LogIO::DEBUG1
718 : << "hertz_p = " << hertz_p.getName()
719 0 : << "\nvalues[0].unit() = " << values[0].unit().getName()
720 0 : << LogIO::POST;
721 0 : }
722 :
723 0 : void FluxCalc_SS_JPL_Butler::compute_venus(Vector<Flux<Double> >& values,
724 : Vector<Flux<Double> >& errors,
725 : const Double angdiam,
726 : const Vector<MFrequency>& mfreqs)
727 : {
728 0 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_venus"));
729 0 : const uInt nfreqs = mfreqs.nelements();
730 0 : Vector<Double> temps(nfreqs);
731 0 : Vector<Float> ghzfreqs(nfreqs);
732 :
733 0 : Float minfreq = 0.303;
734 0 : Float maxfreq = 350.0;
735 0 : for(uInt f = 0; f < nfreqs; ++f){
736 0 : Float ghzfreq = 1.0e-9 * mfreqs[f].get(hertz_p).getValue();
737 :
738 0 : if(ghzfreq < minfreq)
739 0 : minfreq = ghzfreq;
740 0 : else if(ghzfreq > maxfreq)
741 0 : maxfreq = ghzfreq;
742 0 : ghzfreqs[f] = ghzfreq;
743 : }
744 :
745 0 : const uInt nmeas = 75;
746 :
747 : // GHz. Nominally const, but there is no constructor for making a Block from
748 : // a const C array.
749 0 : Float measfreqarr[nmeas] = {350.000, 318.182, 289.256, 262.960, 239.055,
750 : 217.322, 197.566, 179.605, 163.278, 148.434,
751 : 134.940, 122.673, 111.521, 101.383, 92.166,
752 : 83.787, 76.170, 69.246, 62.951, 57.228,
753 : 52.025, 47.296, 42.996, 39.087, 35.534,
754 : 32.304, 29.367, 26.697, 24.270, 22.064,
755 : 20.058, 18.235, 16.577, 15.070, 13.700,
756 : 12.454, 11.322, 10.293, 9.357, 8.507,
757 : 7.733, 7.030, 6.391, 5.810, 5.282,
758 : 4.802, 4.365, 3.968, 3.608, 3.280,
759 : 2.981, 2.710, 2.464, 2.240, 2.036,
760 : 1.851, 1.683, 1.530, 1.391, 1.264,
761 : 1.149, 1.045, 0.950, 0.864, 0.785,
762 : 0.714, 0.649, 0.590, 0.536, 0.487,
763 : 0.443, 0.403, 0.366, 0.333, 0.303};
764 0 : Float *measfreqptr = measfreqarr; // Need a referenceable pointer.
765 0 : const Block<Float> measfreqblk(nmeas, measfreqptr, false);
766 :
767 : // Double to match the type of temps. Nominally const, but there is no
768 : // constructor for making a Block from a const C array.
769 0 : Double meastbarr[nmeas] = {270.2, 273.8, 277.7, 282.0, 286.8,
770 : 292.1, 297.6, 303.5, 309.7, 316.1,
771 : 322.8, 329.6, 336.6, 343.7, 351.1,
772 : 358.7, 366.7, 374.9, 383.6, 392.7,
773 : 402.3, 412.3, 422.8, 433.8, 445.3,
774 : 457.3, 469.8, 482.8, 496.1, 509.9,
775 : 524.1, 538.8, 553.7, 568.7, 583.8,
776 : 598.7, 613.0, 626.5, 638.9, 648.0,
777 : 657.0, 665.0, 671.0, 676.0, 679.0,
778 : 680.0, 680.0, 679.0, 676.0, 671.0,
779 : 664.0, 655.0, 646.0, 638.0, 631.0,
780 : 624.0, 617.0, 610.0, 604.0, 598.0,
781 : 592.0, 586.0, 580.0, 575.0, 570.0,
782 : 565.0, 560.0, 555.0, 550.0, 545.0,
783 : 540.0, 535.0, 530.0, 525.0, 520.0};
784 0 : Double *meastbptr = meastbarr;
785 0 : const Block<Double> meastbblk(nmeas, meastbptr, false);
786 :
787 : // Just let it extrapolate if necessary; the temperature_p given in the
788 : // ephemeris (735K) is so high that I think it's for the surface.
789 0 : InterpolateArray1D<Float, Double>::interpolate(temps, ghzfreqs,
790 0 : Vector<Float>(measfreqblk.begin(),measfreqblk.end()),
791 0 : Vector<Double>(meastbblk.begin(),meastbblk.end()),
792 : InterpolateArray1D<Float, Double>::cubic);
793 :
794 0 : if(minfreq < 0.303)
795 : os << LogIO::WARN
796 : << "At least one of the frequencies, " << minfreq
797 : << "GHz, is below the lower limit for the model (0.303GHz).\n"
798 0 : << LogIO::POST;
799 0 : if(maxfreq > 350.0)
800 : os << LogIO::WARN
801 : << "At least one of the frequencies, " << maxfreq
802 : << "GHz, is above the upper limit for the model (350.0GHz).\n"
803 0 : << LogIO::POST;
804 0 : if(minfreq < 0.303 || maxfreq > 350.0)
805 : os << LogIO::WARN
806 : << "The extrapolation may be very bad."
807 0 : << LogIO::POST;
808 :
809 0 : compute_GB(values, errors, angdiam, mfreqs, temps);
810 0 : }
811 :
812 0 : void FluxCalc_SS_JPL_Butler::compute_jupiter(Vector<Flux<Double> >& values,
813 : Vector<Flux<Double> >& errors,
814 : const Double angdiam,
815 : const Vector<MFrequency>& mfreqs)
816 : {
817 0 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_jupiter"));
818 0 : Bool outOfFreqRange = false;
819 0 : const uInt nfreqs = mfreqs.nelements();
820 0 : Vector<Double> temps(nfreqs);
821 :
822 0 : for(uInt f = 0; f < nfreqs; ++f){
823 0 : Double freq = mfreqs[f].get(hertz_p).getValue();
824 0 : Double lambdacm = 100.0 * C::c / freq; // Wavelength in cm.
825 :
826 0 : if(lambdacm < 0.1){
827 0 : outOfFreqRange = true;
828 0 : lambdacm = 0.1;
829 : }
830 0 : else if(lambdacm > 6.2){
831 0 : outOfFreqRange = true;
832 0 : lambdacm = 6.2;
833 : }
834 :
835 0 : if(lambdacm < 0.44){
836 0 : temps[f] = 170.0;
837 : }
838 0 : else if(lambdacm < 0.7){
839 : // 21.537539 = 10.0 / ln(0.7 / 0.44)
840 0 : temps[f] = 160.0 + 21.537539 * log(0.7 / lambdacm);
841 : }
842 0 : else if(lambdacm < 1.3){
843 : // 48.462196889 = 30.0 / ln(1.3 / 0.7)
844 0 : temps[f] = 130.0 + 48.462196889 * log(1.3 / lambdacm);
845 : }
846 : else
847 : // 65.38532335444 = 100.0 / ln(6.0 / 1.3)
848 0 : temps[f] = 130.0 + 65.38532335444 * log10(lambdacm / 1.3);
849 : }
850 :
851 0 : if(outOfFreqRange)
852 : os << LogIO::WARN
853 : << "At least one of the wavelengths went outside the nominal range\n"
854 : << "of 1mm to 6.2cm, so the wavelength was clamped to 1mm or 6.2cm for\n"
855 : << "calculating the effective temperature of Jupiter."
856 0 : << LogIO::POST;
857 :
858 0 : compute_GB(values, errors, angdiam, mfreqs, temps);
859 0 : }
860 :
861 0 : void FluxCalc_SS_JPL_Butler::compute_uranus(Vector<Flux<Double> >& values,
862 : Vector<Flux<Double> >& errors,
863 : const Double angdiam,
864 : const Vector<MFrequency>& mfreqs)
865 : {
866 0 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_uranus"));
867 0 : Bool outOfFreqRange = false;
868 0 : const uInt nfreqs = mfreqs.nelements();
869 0 : Vector<Double> temps(nfreqs);
870 :
871 0 : for(uInt f = 0; f < nfreqs; ++f){
872 0 : Double freq = mfreqs[f].get(hertz_p).getValue();
873 0 : Double lambdacm = 100.0 * C::c / freq; // Wavelength in cm.
874 :
875 0 : if(lambdacm < 0.07){
876 0 : outOfFreqRange = true;
877 0 : lambdacm = 0.07;
878 : }
879 0 : else if(lambdacm > 6.2){
880 0 : outOfFreqRange = true;
881 0 : lambdacm = 6.2;
882 : }
883 :
884 0 : if(lambdacm < 0.4){
885 : // 32.46063842 = 40.0 / ln(4.0)
886 0 : temps[f] = 90.0 + 32.46063842 * log(10.0 * lambdacm);
887 : }
888 0 : else if(lambdacm < 1.0){
889 0 : temps[f] = 135.0;
890 : }
891 : else
892 0 : temps[f] = 135.0 + 105.0 * log10(lambdacm);
893 : }
894 :
895 0 : if(outOfFreqRange)
896 : os << LogIO::WARN
897 : << "At least one of the wavelengths went outside the nominal range\n"
898 : << "of 0.7mm to 6.2cm, so the wavelength was clamped at either 0.7mm or 6.2cm\n"
899 : << "for calculating the effective temperature of Uranus."
900 0 : << LogIO::POST;
901 :
902 0 : compute_GB(values, errors, angdiam, mfreqs, temps);
903 0 : }
904 :
905 0 : void FluxCalc_SS_JPL_Butler::compute_neptune(Vector<Flux<Double> >& values,
906 : Vector<Flux<Double> >& errors,
907 : const Double angdiam,
908 : const Vector<MFrequency>& mfreqs)
909 : {
910 0 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_neptune"));
911 0 : Bool outOfFreqRange = false;
912 0 : const uInt nfreqs = mfreqs.nelements();
913 0 : Vector<Double> temps(nfreqs);
914 :
915 0 : for(uInt f = 0; f < nfreqs; ++f){
916 0 : Double freq = 1.0e-9 * mfreqs[f].get(hertz_p).getValue(); // GHz
917 :
918 : // These temperatures agree well with the ones collected in
919 : // http://adsabs.harvard.edu/abs/1995EM&P...67...89S (Spilker)
920 : // (0.1-20cm = 1.5-300GHz) and
921 : // http://www.ericweisstein.com/research/papers/applopt/node10.html
922 : // (200-300GHz) except around 70 GHz. It's not so much a knee
923 : // at 70 GHz as a flat interval around 70 GHz. Unfortunately I don't yet
924 : // have any data right there, but Planck should publish better models for
925 : // the planets.
926 : //
927 : // (Spilker attributed the flat interval or dip in Jupiter, Saturn, and
928 : // Neptune to NH3. Uranus has barely any, at least while we're looking at
929 : // its pole.)
930 : //
931 0 : if(freq < 4.0){
932 0 : outOfFreqRange = true;
933 0 : freq = 4.0;
934 : }
935 0 : else if(freq > 1000.0){
936 0 : outOfFreqRange = true;
937 0 : freq = 1000.0;
938 : }
939 :
940 0 : if(freq < 70.0){
941 : // 30.083556662 = 80.0 / ln(1000.0 / 70.0)
942 0 : temps[f] = 140.0 - 30.083556662 * log(freq / 70.0);
943 : }
944 : else
945 : // 34.93815 = 100.0 / ln(70.0 / 4.0)
946 0 : temps[f] = 240.0 - 34.93815 * log(freq / 4.0);
947 : }
948 :
949 0 : if(outOfFreqRange)
950 : os << LogIO::WARN
951 : << "At least one of the frequencies went outside the nominal range\n"
952 : << "of 4 to 1000 GHz for Neptune, so it was clamped to 4 or 1000 GHz\n"
953 : << "for calculating the effective temperature."
954 0 : << LogIO::POST;
955 :
956 0 : compute_GB(values, errors, angdiam, mfreqs, temps);
957 0 : }
958 :
959 14 : Bool FluxCalc_SS_JPL_Butler::compute_constant_temperature(Vector<Flux<Double> >& values,
960 : Vector<Flux<Double> >& errors,
961 : const Double angdiam,
962 : const Vector<MFrequency>& mfreqs,
963 : const Bool report)
964 : {
965 28 : LogIO os(LogOrigin("FluxCalc_SS_JPL_Butler", "compute_constant_temperature"));
966 :
967 14 : Bool success = true;
968 14 : Double ephem_temp = temperature_p; // Store it.
969 :
970 14 : switch(objnum_p){
971 0 : case FluxCalc_SS_JPL_Butler::Pluto:
972 0 : if(report)
973 : os << LogIO::NORMAL1
974 : << "Using the value from:"
975 : << " Altenhoff, W. J., R. Chini, H. Hein, E. Kreysa, P. G. Mezger, "
976 : << " C. Salter, and J. B. Schraml, First radio astronomical estimate "
977 : << " of the temperature of Pluto, A&ALett, 190, L15-L17, 1988"
978 : << "which is: Tb = 35 K at 1.27 mm. this is a correction from the "
979 : << "value of 39 K in the paper, due to the incorrect geometric mean size "
980 : << "used for Pluto and Charon (1244 km vs the correct 1320 km). this is"
981 : << "similar to the value found in:"
982 : << " Stern, S. A., D. A. Weintraub, and M. C. Festou, Evidence for a Low "
983 : << " Surface Temperature on Pluto from Millimeter-Wave Thermal"
984 : << " Emission Measurements, Science, 261, 1713-1716, 1993"
985 : << "and is a good match to the physical temperature reported in:"
986 : << " Tryka, K. A., R. H. Brown, D. P. Cruikshank, T. C. Owen, T. R."
987 : << " Geballe, and C. DeBergh, Temperature of Nitrogen Ice on Pluto and"
988 : << " Its Implications for Flux Measurements, Icarus, 112, 513-527, "
989 : << " 1994"
990 : << "where they give a surface temperature of 40+-2 K. this would imply"
991 : << "an emissivity of 0.875, which is certainly reasonable."
992 0 : << LogIO::POST;
993 0 : temperature_p = 35.0;
994 0 : break;
995 0 : case FluxCalc_SS_JPL_Butler::Io:
996 0 : if(report)
997 : os << LogIO::NORMAL1
998 : << "Reference for Io's temperature (110K):\n"
999 : << "Rathbun, J.A., Spencer, J.R. Tamppari, L.K., Martin, T.Z., Barnard, L.,\n"
1000 : << "Travis, L.D. (2004). \"Mapping of Io's thermal radiation by the Galileo\n"
1001 : << "photopolarimeter-radiometer (PPR) instrument\". Icarus 169:127-139.\n"
1002 : << "doi:10.1016/j.icarus.2003.12.021\n"
1003 0 : << LogIO::POST;
1004 0 : temperature_p = 110.0;
1005 0 : break;
1006 0 : case FluxCalc_SS_JPL_Butler::Ganymede:
1007 0 : if(report)
1008 : os << LogIO::NORMAL1
1009 : << "Reference for Ganymede's temperature (110K):\n"
1010 : << "Delitsky, Mona L., Lane, Arthur L. (1998). \"Ice chemistry of Galilean satellites\"\n"
1011 : << "J. of Geophys. Res. 103 (E13): 31,391-31,403. doi:10.1029/1998/JE900020\n"
1012 : << "http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/20675/1/98-1725.pdf"
1013 0 : << LogIO::POST;
1014 0 : temperature_p = 110.0;
1015 0 : break;
1016 0 : case FluxCalc_SS_JPL_Butler::Callisto:
1017 0 : if(report)
1018 : os << LogIO::NORMAL1
1019 : << "Reference for Callisto's temperature (134 +- 11 K):\n"
1020 : << "Moore, Jeffrey M., Chapman, Clark R., Bierhaus, Edward B. et al\n"
1021 : << "(2004). \"Callisto\" In Bagenal, F., Dowling, T.E., McKinnon, W.B.,\n"
1022 : << "\"Jupiter: The Planet, Satellites, and Magnetosphere\". Cambridge Univ. Press"
1023 0 : << LogIO::POST;
1024 0 : temperature_p = 134.0; // +- 11.
1025 0 : break;
1026 0 : case FluxCalc_SS_JPL_Butler::Europa:
1027 0 : if(report)
1028 : os << LogIO::NORMAL1
1029 : << "Reference for Europa's temperature (109 K):\n"
1030 : << "http://science.nasa.gov/science-news/science-at-nasa/1998/ast03dec98_1/"
1031 0 : << LogIO::POST;
1032 0 : temperature_p = 109.0;
1033 0 : break;
1034 14 : case FluxCalc_SS_JPL_Butler::Titan:
1035 14 : temperature_p = 76.6;
1036 14 : break;
1037 0 : case FluxCalc_SS_JPL_Butler::Triton:
1038 0 : if(report)
1039 : os << LogIO::NORMAL1
1040 : << "Reference for Triton's temperature (38 K):\n"
1041 : << "http://solarsystem.nasa.gov/planets/profile.cfm?Object=Triton"
1042 0 : << LogIO::POST;
1043 0 : temperature_p = 38.0;
1044 0 : break;
1045 0 : case FluxCalc_SS_JPL_Butler::Ceres:
1046 0 : if(report)
1047 : os << LogIO::NORMAL1
1048 : << "Reference for Ceres' mean temperature (167K):\n"
1049 : << "Saint-Pe, O., Combes, N., Rigaut, F. (1993). \"Ceres surface properties\n"
1050 : << "by high-resolution imaging from Earth\" Icarus 105:271-281.\n"
1051 : << "doi:10.1006/icar.1993.1125"
1052 0 : << LogIO::POST;
1053 0 : temperature_p = 167.0;
1054 0 : break;
1055 0 : case FluxCalc_SS_JPL_Butler::Pallas:
1056 0 : if(report)
1057 : os << LogIO::WARN
1058 : << "The orbit of Pallas has an eccentricity of 0.231. This is not yet accounted\n"
1059 : << "for when calculating its temperature (taken to be 164K)."
1060 0 : << LogIO::POST;
1061 0 : temperature_p = 164.0;
1062 0 : break;
1063 0 : case FluxCalc_SS_JPL_Butler::Juno:
1064 0 : if(report){
1065 : os << LogIO::WARN
1066 : << "Juno has a large crater and temperature changes that CASA does not fully account for."
1067 0 : << LogIO::POST;
1068 : os << LogIO::NORMAL1
1069 : << "Reference for Juno's mean temperature (163K):\n"
1070 : << "Lim, Lucy F., McConnochie, Timothy H., Bell, James F., Hayward, Thomas L. (2005).\n"
1071 : << "\"Thermal infrared (8-13 um) spectra of 29 asteroids: the Cornell Mid-Infrared\n"
1072 : << "Asteroid Spectroscopy (MIDAS) Survey\" Icarus 173:385-408.\n"
1073 : << "doi:10.1016/j.icarus.2004.08.005."
1074 0 : << LogIO::POST;
1075 : }
1076 0 : temperature_p = 163.0;
1077 0 : break;
1078 0 : case FluxCalc_SS_JPL_Butler::Vesta:
1079 0 : if(report){
1080 : os << LogIO::WARN
1081 : << "Vesta has a large crater, and its mean emissivity varies from\n"
1082 : << "0.6 in the submm to 0.7 in the mm. Its mean submm brightness temperature varies\n"
1083 : << "from ~116 to 173K and its rotation period is 5.342h.\n "
1084 : << "CASA does not yet account for these variations."
1085 0 : << LogIO::POST;
1086 : os << LogIO::NORMAL1
1087 : << "Reference for Vesta's mean brightness temperature (160K):\n"
1088 : << "Chamberlain et al., (2009).\n"
1089 : << "\"Submillimeter photometry and lightcurves of Ceres and other large asteroids\"\n"
1090 : << "Icarus 202:487-501.\n"
1091 0 : << LogIO::POST;
1092 : }
1093 0 : temperature_p = 160.0;
1094 0 : break;
1095 0 : default:
1096 0 : break;
1097 : };
1098 :
1099 14 : if(temperature_p > 0.0){
1100 14 : if(report)
1101 4 : os << LogIO::NORMAL << "Using blackbody model." << LogIO::POST;
1102 14 : compute_BB(values, errors, angdiam, mfreqs);
1103 : }
1104 : else{
1105 : os << LogIO::SEVERE
1106 : << "An ephemeris was found, but not a temperature."
1107 0 : << LogIO::POST;
1108 0 : success = false;
1109 : }
1110 :
1111 14 : if(ephem_temp > 0.0)
1112 10 : temperature_p = ephem_temp; // Restore it.
1113 :
1114 28 : return success;
1115 : }
1116 :
1117 0 : Bool FluxCalc_SS_JPL_Butler::setObj(const String& objname)
1118 : {
1119 0 : name_p = objname;
1120 0 : return setObjNum();
1121 : }
1122 :
1123 : } //# NAMESPACE CASA - END
1124 :
|
