00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include <iostream>
00036 #include <string>
00037 #include <sstream>
00038 #include <vector>
00039 #include <map>
00040 #include <cstdlib>
00041 #include <cctype>
00042 #include <cstring>
00043 #include "fitsio.h"
00044 #include "fitshandle.h"
00045 #include "cxxutils.h"
00046
00047 using namespace std;
00048
00049 namespace {
00050
00051 string ftc2char (int type)
00052 {
00053 switch (type)
00054 {
00055 case TLOGICAL : return "L";
00056 case TFLOAT : return "E";
00057 case TDOUBLE : return "D";
00058 case TBYTE : return "B";
00059 case TSHORT : return "I";
00060 case TINT32BIT: return "J";
00061 case TLONGLONG: return "K";
00062 case TSTRING : return "A";
00063 default: throw Message_error("wrong datatype in ftc2char()");
00064 }
00065 }
00066
00067 string ftc2asciiform (int type)
00068 {
00069 switch (type)
00070 {
00071 case TFLOAT: return "E14.7";
00072 case TDOUBLE: return "D23.15";
00073 case TBYTE: return "I4";
00074 case TSHORT: return "I6";
00075 case TINT32BIT: return "I11";
00076 case TLONGLONG: return "I22";
00077 default: throw Message_error("wrong datatype in ftc2asciiform()");
00078 }
00079 }
00080
00081 string fixkey (const string &key)
00082 {
00083 for (unsigned int m=0; m<key.size(); ++m)
00084 if (islower(key[m])) return string("HIERARCH "+key);
00085
00086 return key;
00087 }
00088 }
00089
00090 void fitshandle::check_errors() const
00091 {
00092 if (status==0) return;
00093 char msg[81];
00094 fits_get_errstatus (status, msg);
00095 cerr << msg << endl;
00096 while (fits_read_errmsg(msg)) cerr << msg << endl;
00097 throw Message_error("FITS error");
00098 }
00099
00100 void fitshandle::clean_data()
00101 {
00102 if (!fptr) return;
00103 axes_.clear();
00104 columns_.clear();
00105 hdutype_=INVALID;
00106 bitpix_=INVALID;
00107 nrows_=0;
00108 }
00109
00110 void fitshandle::clean_all()
00111 {
00112 if (!fptr) return;
00113 clean_data();
00114 fits_close_file (fptr, &status);
00115 check_errors();
00116 fptr=0;
00117 }
00118
00119 void fitshandle::init_image()
00120 {
00121 int naxis;
00122 fits_get_img_type (fptr, &bitpix_, &status);
00123 fits_get_img_dim (fptr, &naxis, &status);
00124 check_errors();
00125 arr<LONGLONG> naxes(naxis);
00126 fits_get_img_sizell (fptr, naxis, &naxes[0], &status);
00127 for (long m=0; m<naxis; ++m) axes_.push_back(naxes[naxis-m-1]);
00128 check_errors();
00129 }
00130
00131 void fitshandle::init_asciitab()
00132 {
00133 char ttype[81], tunit[81], tform[81];
00134 int ncol, typecode;
00135 fits_get_num_cols (fptr, &ncol, &status);
00136 { LONGLONG tmp; fits_get_num_rowsll (fptr, &tmp, &status); nrows_=tmp; }
00137 check_errors();
00138 for (int m=1; m<=ncol; ++m)
00139 {
00140 fits_get_acolparms (fptr, m, ttype, 0, tunit, tform,
00141 0, 0, 0, 0, &status);
00142 fits_ascii_tform (tform, &typecode, 0,0, &status);
00143 check_errors();
00144 columns_.push_back (fitscolumn (ttype,tunit,1,typecode));
00145 }
00146 }
00147
00148 void fitshandle::init_bintab()
00149 {
00150 char ttype[81], tunit[81], tform[81];
00151 LONGLONG repc;
00152 int ncol, typecode;
00153 fits_get_num_cols (fptr, &ncol, &status);
00154 { LONGLONG tmp; fits_get_num_rowsll (fptr, &tmp, &status); nrows_=tmp; }
00155 check_errors();
00156 for (int m=1; m<=ncol; ++m)
00157 {
00158 fits_get_bcolparmsll (fptr, m, ttype, tunit, tform, &repc,
00159 0, 0, 0, 0, &status);
00160 fits_binary_tform (tform, &typecode, 0,0, &status);
00161 check_errors();
00162 columns_.push_back (fitscolumn (ttype,tunit,repc,typecode));
00163 }
00164 }
00165
00166 void fitshandle::init_data()
00167 {
00168 clean_data();
00169 fits_get_hdu_type (fptr, &hdutype_, &status);
00170 check_errors();
00171 switch(hdutype_)
00172 {
00173 case IMAGE_HDU:
00174 init_image();
00175 break;
00176 case ASCII_TBL:
00177 init_asciitab();
00178 break;
00179 case BINARY_TBL:
00180 init_bintab();
00181 break;
00182 default:
00183 throw Message_error("init_data(): wrong HDU type");
00184 break;
00185 }
00186 }
00187
00188 void fitshandle::read_col (int colnum, void *data, int64 ndata, int dtype,
00189 int64 offset) const
00190 {
00191 assert_table_hdu("fitshandle::read_column()",colnum);
00192 int64 repc = columns_[colnum-1].repcount();
00193 planck_assert (ndata<=(repc*nrows_-offset),
00194 "read_column(): array too large");
00195 int64 frow = offset/repc+1;
00196 int64 felem = offset%repc+1;
00197 fits_read_col (fptr, dtype, colnum, frow, felem, ndata, 0, data, 0, &status);
00198 check_errors();
00199 }
00200 void fitshandle::write_col (int colnum, const void *data, int64 ndata,
00201 int dtype, int64 offset)
00202 {
00203 assert_table_hdu("fitshandle::write_column()",colnum);
00204 int64 repc = columns_[colnum-1].repcount();
00205 int64 frow = offset/repc+1;
00206 int64 felem = offset%repc+1;
00207 fits_write_col (fptr, dtype, colnum, frow, felem, ndata,
00208 const_cast<void *>(data), &status);
00209 nrows_ = max(nrows_,offset+ndata);
00210 check_errors();
00211 }
00212
00213 void fitshandle::open (const string &fname, int rwmode)
00214 {
00215 clean_all();
00216 fits_open_file(&fptr, fname.c_str(), rwmode, &status);
00217 check_errors();
00218 init_data();
00219 }
00220
00221 void fitshandle::create (const string &fname)
00222 {
00223 clean_all();
00224 fits_create_file(&fptr, fname.c_str(), &status);
00225 fits_write_imghdr(fptr, 8, 0, 0, &status);
00226 fits_write_date(fptr, &status);
00227 check_errors();
00228 init_data();
00229 }
00230
00231
00232 void fitshandle::delete_file (const string &name)
00233 {
00234 fitsfile *ptr;
00235 int status = 0;
00236 fits_open_file(&ptr, name.c_str(), READWRITE, &status);
00237 fits_delete_file(ptr, &status);
00238 if (status==0) return;
00239
00240 char msg[81];
00241 fits_get_errstatus (status, msg);
00242 cerr << msg << endl;
00243 while (fits_read_errmsg(msg)) cerr << msg << endl;
00244 throw Message_error("FITS error");
00245 }
00246
00247 void fitshandle::goto_hdu (int hdu)
00248 {
00249 int curhdu;
00250 fits_get_hdu_num(fptr,&curhdu);
00251 if (curhdu!=hdu)
00252 {
00253 fits_movabs_hdu(fptr, hdu, &hdutype_, &status);
00254 check_errors();
00255 init_data();
00256 }
00257 }
00258
00259 int fitshandle::num_hdus () const
00260 {
00261 int result;
00262 fits_get_num_hdus (fptr, &result, &status);
00263 check_errors();
00264 return result;
00265 }
00266
00267 void fitshandle::insert_bintab (const vector<fitscolumn> &cols,
00268 const string &extname)
00269 {
00270 clean_data();
00271 int ncol=cols.size();
00272 arr2b<char> ttype(ncol,81), tform(ncol,81), tunit(ncol,81);
00273
00274 for (long m=0; m<ncol; ++m)
00275 {
00276 strcpy (ttype[m], cols[m].name().c_str());
00277 strcpy (tunit[m], cols[m].unit().c_str());
00278 ostringstream x;
00279 x << cols[m].repcount() << ftc2char (cols[m].type());
00280 strcpy (tform[m], x.str().c_str());
00281 }
00282 fits_insert_btbl (fptr, nrows_, ncol, ttype.p0(), tform.p0(), tunit.p0(),
00283 const_cast<char *>(extname.c_str()), 0, &status);
00284 check_errors();
00285 init_data();
00286 }
00287
00288 void fitshandle::insert_asctab (const vector<fitscolumn> &cols,
00289 const string &extname)
00290 {
00291 clean_data();
00292 int ncol=cols.size();
00293 arr2b<char> ttype(ncol,81), tform(ncol,81), tunit(ncol,81);
00294
00295 for (long m=0; m<ncol; ++m)
00296 {
00297 strcpy (ttype[m], cols[m].name().c_str());
00298 strcpy (tunit[m], cols[m].unit().c_str());
00299 ostringstream x;
00300 if (cols[m].type()!=TSTRING)
00301 {
00302 planck_assert (cols[m].repcount()==1,"bad repcount for ASCII table");
00303 x << ftc2asciiform (cols[m].type());
00304 }
00305 else
00306 {
00307 x << "A" << dataToString(cols[m].repcount());
00308 }
00309 strcpy (tform[m], x.str().c_str());
00310 }
00311 fits_insert_atbl (fptr, 0, nrows_, ncol, ttype.p0(), 0, tform.p0(),
00312 tunit.p0(), const_cast<char *>(extname.c_str()), &status);
00313 check_errors();
00314 init_data();
00315 }
00316
00317 void fitshandle::insert_image (int btpx, const vector<int64> &Axes)
00318 {
00319 clean_data();
00320 arr<LONGLONG> tmpax(Axes.size());
00321 for (long m=0; m<long(Axes.size()); m++) tmpax[m]=Axes[Axes.size()-1-m];
00322 fits_insert_imgll (fptr, btpx, Axes.size(), &tmpax[0], &status);
00323 check_errors();
00324 init_data();
00325 }
00326
00327 template<typename T>
00328 void fitshandle::insert_image (int btpx, const arr2<T> &data)
00329 {
00330 clean_data();
00331 arr<LONGLONG> tmpax(2);
00332 tmpax[0] = data.size2(); tmpax[1] = data.size1();
00333 fits_insert_imgll (fptr, btpx, 2, &tmpax[0], &status);
00334 arr2<T> &tmparr = const_cast<arr2<T> &> (data);
00335 fits_write_img (fptr, FITSUTIL<T>::DTYPE, 1, tmpax[0]*tmpax[1],
00336 &tmparr[0][0], &status);
00337 check_errors();
00338 init_data();
00339 }
00340
00341 template void fitshandle::insert_image (int btpx, const arr2<double> &data);
00342 template void fitshandle::insert_image (int btpx, const arr2<float> &data);
00343 template void fitshandle::insert_image (int btpx, const arr2<int> &data);
00344
00345 void fitshandle::write_checksum()
00346 {
00347 assert_connected("fitshandle::write_checksum()");
00348 fits_write_chksum (fptr, &status);
00349 check_errors();
00350 }
00351
00352 void fitshandle::copy_historified_header (const fitshandle &orig)
00353 {
00354 const char *inclist[] = { "*" };
00355 const char *exclist[] = {
00356 "SIMPLE","BITPIX","NAXIS","NAXIS#","PCOUNT","GCOUNT",
00357 "EXTEND","ORIGIN","DATE*","TFIELDS","TTYPE#","TFORM#",
00358 "TUNIT#","EXTNAME","CTYPE#","CRVAL#","CRPIX#","CDELT#",
00359 "XTENSION","INSTRUME","TELESCOP","PDMTYPE","TBCOL#" };
00360 char card[81];
00361 string card2;
00362 orig.assert_connected("fitshandle::copy_historified_header()");
00363 assert_connected("fitshandle::copy_historified_header()");
00364 fits_read_record (orig.fptr, 0, card, &status);
00365 check_errors();
00366 while (true)
00367 {
00368 fits_find_nextkey (orig.fptr, const_cast<char **>(inclist), 1,
00369 const_cast<char **>(exclist), 23, card, &status);
00370 if (status!=0) break;
00371 card2=trim(card);
00372 if (card2!="END" && card2!="COMMENT" && card2!="HISTORY")
00373 {
00374 if (card2.find("COMMENT")==0)
00375 card2.replace(0,7,"HISTORY");
00376 if (card2.find("HISTORY")!=0)
00377 card2.insert(0,"HISTORY ");
00378 if (card2.length()<=80)
00379 fits_write_record (fptr, card2.c_str(), &status);
00380 else
00381 {
00382 fits_write_record (fptr, card2.substr(0,80).c_str(), &status);
00383 card2=card2.substr(80,string::npos);
00384 card2.insert(0,"HISTORY ");
00385 fits_write_record (fptr, card2.c_str(),&status);
00386 }
00387 }
00388 check_errors();
00389 }
00390 if (status==KEY_NO_EXIST) { fits_clear_errmsg(); status=0; }
00391 check_errors();
00392 }
00393
00394 void fitshandle::copy_header (const fitshandle &orig)
00395 {
00396 const char *inclist[] = { "*" };
00397 const char *exclist[] = {
00398 "SIMPLE","BITPIX","NAXIS","NAXIS#","PCOUNT","GCOUNT",
00399 "EXTEND","ORIGIN","DATE*","TFIELDS","TTYPE#","TFORM#",
00400 "TUNIT#","EXTNAME","CTYPE#","CRVAL#","CRPIX#","CDELT#",
00401 "XTENSION","INSTRUME","TELESCOP","PDMTYPE","TBCOL#" };
00402 char card[81];
00403 string card2;
00404 orig.assert_connected("fitshandle::copy_header()");
00405 assert_connected("fitshandle::copy_header()");
00406 fits_read_record (orig.fptr, 0, card, &status);
00407 check_errors();
00408 while (true)
00409 {
00410 fits_find_nextkey (orig.fptr, const_cast<char **>(inclist), 1,
00411 const_cast<char **>(exclist), 23, card, &status);
00412 if (status!=0) break;
00413 card2=trim(card);
00414 if (card2!="END" && card2!="COMMENT" && card2!="HISTORY")
00415 fits_write_record (fptr, card, &status);
00416 check_errors();
00417 }
00418 if (status==KEY_NO_EXIST) { fits_clear_errmsg(); status=0; }
00419 check_errors();
00420 }
00421
00422 void fitshandle::get_all_keys(vector<string> &keys) const
00423 {
00424 keys.clear();
00425 char card[81];
00426 const char *inclist[] = { "*" };
00427 assert_connected("fitshandle::get_all_keys()");
00428 fits_read_record (fptr, 0, card, &status);
00429 check_errors();
00430 while (true)
00431 {
00432 fits_find_nextkey (fptr, const_cast<char **>(inclist), 1,
00433 0, 0, card, &status);
00434 if (status!=0) break;
00435 if (fits_get_keyclass(card)==TYP_USER_KEY)
00436 {
00437 char keyname[80];
00438 int dummy;
00439 fits_get_keyname(card, keyname, &dummy, &status);
00440 check_errors();
00441 keys.push_back(keyname);
00442 }
00443 check_errors();
00444 }
00445 if (status==KEY_NO_EXIST) { fits_clear_errmsg(); status=0; }
00446 check_errors();
00447 }
00448
00449 void fitshandle::check_key_present(const string &name) const
00450 {
00451 char card[81];
00452 fits_read_card(fptr, const_cast<char *>(name.c_str()), card, &status);
00453 if (status==KEY_NO_EXIST)
00454 { fits_clear_errmsg(); status=0; return; }
00455 check_errors();
00456
00457
00458 }
00459
00460 template<typename T> void fitshandle::add_key (const string &name,
00461 const T &value, const string &comment)
00462 {
00463 assert_connected("fitshandle::add_key()");
00464 string name2 = fixkey(name);
00465 check_key_present (name);
00466 fits_write_key (fptr, FITSUTIL<T>::DTYPE, const_cast<char *>(name2.c_str()),
00467 const_cast<T *>(&value), const_cast<char *>(comment.c_str()), &status);
00468 check_errors();
00469 }
00470
00471 template void fitshandle::add_key(const string &name,
00472 const signed char &value, const string &comment);
00473 template void fitshandle::add_key(const string &name,
00474 const short &value, const string &comment);
00475 template void fitshandle::add_key(const string &name,
00476 const int &value, const string &comment);
00477 template void fitshandle::add_key(const string &name,
00478 const long &value, const string &comment);
00479 template void fitshandle::add_key(const string &name,
00480 const long long &value, const string &comment);
00481 template void fitshandle::add_key(const string &name,
00482 const float &value, const string &comment);
00483 template void fitshandle::add_key(const string &name,
00484 const double &value, const string &comment);
00485 template<> void fitshandle::add_key(const string &name,
00486 const bool &value, const string &comment)
00487 {
00488 assert_connected("fitshandle::add_key()");
00489 string name2 = fixkey(name);
00490 check_key_present (name);
00491 int val=value;
00492 fits_write_key (fptr, TLOGICAL, const_cast<char *>(name2.c_str()),
00493 &val, const_cast<char *>(comment.c_str()),
00494 &status);
00495 check_errors();
00496 }
00497 template<> void fitshandle::add_key (const string &name,
00498 const string &value, const string &comment)
00499 {
00500 assert_connected("fitshandle::add_key()");
00501 string name2 = fixkey(name);
00502 check_key_present (name);
00503 fits_write_key_longstr (fptr, const_cast<char *>(name2.c_str()),
00504 const_cast<char *>(value.c_str()), const_cast<char *>(comment.c_str()),
00505 &status);
00506 check_errors();
00507 }
00508
00509 template<typename T> void fitshandle::update_key (const string &name,
00510 const T &value, const string &comment)
00511 {
00512 assert_connected("fitshandle::update_key()");
00513 string name2 = fixkey(name);
00514 fits_update_key (fptr, FITSUTIL<T>::DTYPE, const_cast<char *>(name2.c_str()),
00515 const_cast<T *>(&value), const_cast<char *>(comment.c_str()), &status);
00516 check_errors();
00517 }
00518
00519 template void fitshandle::update_key(const string &name,
00520 const signed char &value, const string &comment);
00521 template void fitshandle::update_key(const string &name,
00522 const short &value, const string &comment);
00523 template void fitshandle::update_key(const string &name,
00524 const int &value, const string &comment);
00525 template void fitshandle::update_key(const string &name,
00526 const long &value, const string &comment);
00527 template void fitshandle::update_key(const string &name,
00528 const long long &value, const string &comment);
00529 template void fitshandle::update_key(const string &name,
00530 const float &value, const string &comment);
00531 template void fitshandle::update_key(const string &name,
00532 const double &value, const string &comment);
00533 template<> void fitshandle::update_key(const string &name,
00534 const bool &value, const string &comment)
00535 {
00536 assert_connected("fitshandle::update_key()");
00537 string name2 = fixkey(name);
00538 int val=value;
00539 fits_update_key (fptr, TLOGICAL, const_cast<char *>(name2.c_str()),
00540 &val, const_cast<char *>(comment.c_str()),
00541 &status);
00542 check_errors();
00543 }
00544 template<> void fitshandle::update_key (const string &name,
00545 const string &value, const string &comment)
00546 {
00547 assert_connected("fitshandle::update_key()");
00548 string name2 = fixkey(name);
00549 fits_update_key_longstr (fptr, const_cast<char *>(name2.c_str()),
00550 const_cast<char *>(value.c_str()), const_cast<char *>(comment.c_str()),
00551 &status);
00552 check_errors();
00553 }
00554
00555 void fitshandle::delete_key (const string &name)
00556 {
00557 assert_connected("fitshandle::delete_key()");
00558 fits_delete_key (fptr, const_cast<char *>(name.c_str()), &status);
00559 check_errors();
00560 }
00561
00562 void fitshandle::add_comment(const string &comment)
00563 {
00564 assert_connected("fitshandle::add_comment()");
00565 fits_write_comment(fptr,const_cast<char *>(comment.c_str()),&status);
00566 check_errors();
00567 }
00568
00569 template<typename T> void fitshandle::get_key
00570 (const string &name, T &value) const
00571 {
00572 assert_connected("fitshandle::get_key()");
00573 fits_read_key (fptr, FITSUTIL<T>::DTYPE, const_cast<char *>(name.c_str()),
00574 &value, 0, &status);
00575 if (status==KEY_NO_EXIST) throw Message_error
00576 ("Fitshandle::get_key(): key "+name+" not found");
00577 check_errors();
00578 }
00579 template void fitshandle::get_key(const string &name,signed char &value) const;
00580 template void fitshandle::get_key(const string &name,short &value) const;
00581 template void fitshandle::get_key(const string &name,int &value) const;
00582 template void fitshandle::get_key(const string &name,long &value) const;
00583 template void fitshandle::get_key(const string &name,long long &value) const;
00584 template void fitshandle::get_key(const string &name,float &value) const;
00585 template void fitshandle::get_key(const string &name,double &value) const;
00586 template<> void fitshandle::get_key(const string &name,bool &value) const
00587 {
00588 assert_connected("fitshandle::get_key()");
00589 int val;
00590 fits_read_key (fptr, TLOGICAL, const_cast<char *>(name.c_str()), &val, 0,
00591 &status);
00592 if (status==KEY_NO_EXIST) throw Message_error
00593 ("Fitshandle::get_key(): key "+name+" not found");
00594 check_errors();
00595 value=val;
00596 }
00597 template<> void fitshandle::get_key (const string &name,string &value) const
00598 {
00599 char *tmp=0;
00600 assert_connected("fitshandle::get_key()");
00601 fits_read_key_longstr (fptr, const_cast<char *>(name.c_str()), &tmp, 0,
00602 &status);
00603 if (status==KEY_NO_EXIST) throw Message_error
00604 ("Fitshandle::get_key(): key "+name+" not found");
00605 check_errors();
00606 value=tmp;
00607 if (tmp) free(tmp);
00608 }
00609
00610 bool fitshandle::key_present(const string &name) const
00611 {
00612 char card[81];
00613 assert_connected("fitshandle::key_present()");
00614 fits_read_card(fptr, const_cast<char *>(name.c_str()), card, &status);
00615 if (status==KEY_NO_EXIST)
00616 { fits_clear_errmsg(); status=0; return false; }
00617 check_errors();
00618 return true;
00619 }
00620
00621 int fitshandle::get_key_type(const string &name) const
00622 {
00623 char card[81],value[81],dtype[10];
00624 assert_connected("fitshandle::get_key_type()");
00625 fits_read_card(fptr, const_cast<char *>(name.c_str()), card, &status);
00626 check_errors();
00627 fits_parse_value(card,value,0,&status);
00628 fits_get_keytype(value,dtype,&status);
00629 check_errors();
00630 switch(dtype[0])
00631 {
00632 case 'C' : return PLANCK_STRING;
00633 case 'L' : return PLANCK_BOOL;
00634 case 'I' : return PLANCK_INT64;
00635 case 'F' : return PLANCK_FLOAT64;
00636 default : throw Message_error ("unknown key type");
00637 }
00638 }
00639
00640 void fitshandle::assert_pdmtype (const string &pdmtype) const
00641 {
00642 string type;
00643 get_key("PDMTYPE",type);
00644 if (pdmtype==type) return;
00645 cerr << "PDMTYPE " << pdmtype << " expected, but found " << type << endl;
00646 }
00647
00648 void fitshandle::read_column_raw_void
00649 (int colnum, void *data, int type, int64 num, int64 offset) const
00650 {
00651 switch (type)
00652 {
00653 case PLANCK_INT8:
00654 read_col (colnum, data, num, TBYTE, offset); break;
00655 case PLANCK_INT16:
00656 read_col (colnum, data, num, TSHORT, offset); break;
00657 case PLANCK_INT32:
00658 read_col (colnum, data, num, TINT, offset); break;
00659 case PLANCK_INT64:
00660 read_col (colnum, data, num, TLONGLONG, offset); break;
00661 case PLANCK_FLOAT32:
00662 read_col (colnum, data, num, TFLOAT, offset); break;
00663 case PLANCK_FLOAT64:
00664 read_col (colnum, data, num, TDOUBLE, offset); break;
00665 case PLANCK_BOOL:
00666 read_col (colnum, data, num, TLOGICAL, offset); break;
00667 case PLANCK_STRING:
00668 {
00669 string *data2 = static_cast<string *> (data);
00670 assert_table_hdu("fitshandle::read_column()",colnum);
00671 planck_assert (num<=(nrows_-offset),
00672 "read_column(): array too large");
00673 arr2b<char> tdata(num, columns_[colnum-1].repcount()+1);
00674 fits_read_col (fptr, TSTRING, colnum, offset+1, 1, num,
00675 0, tdata.p0(), 0, &status);
00676 check_errors();
00677 for (long m=0;m<num;++m) data2[m]=tdata[m];
00678 break;
00679 }
00680 default:
00681 throw Message_error ("unsupported data type in read_column_raw_void()");
00682 }
00683 }
00684
00685 void fitshandle::write_column_raw_void
00686 (int colnum, const void *data, int type, int64 num, int64 offset)
00687 {
00688 switch (type)
00689 {
00690 case PLANCK_INT8:
00691 write_col (colnum, data, num, TBYTE, offset); break;
00692 case PLANCK_INT16:
00693 write_col (colnum, data, num, TSHORT, offset); break;
00694 case PLANCK_INT32:
00695 write_col (colnum, data, num, TINT, offset); break;
00696 case PLANCK_INT64:
00697 write_col (colnum, data, num, TLONGLONG, offset); break;
00698 case PLANCK_FLOAT32:
00699 write_col (colnum, data, num, TFLOAT, offset); break;
00700 case PLANCK_FLOAT64:
00701 write_col (colnum, data, num, TDOUBLE, offset); break;
00702 case PLANCK_BOOL:
00703 write_col (colnum, data, num, TLOGICAL, offset); break;
00704 case PLANCK_STRING:
00705 {
00706 const string *data2 = static_cast<const string *> (data);
00707 assert_table_hdu("fitshandle::write_column()",colnum);
00708 int stringlen = columns_[colnum-1].repcount()+1;
00709 arr2b<char> tdata(num, stringlen);
00710 for (long m=0;m<num;++m)
00711 {
00712 strncpy(tdata[m],data2[m].c_str(),stringlen-1);
00713 tdata[m][stringlen-1] = '\0';
00714 }
00715 fits_write_col (fptr, TSTRING, colnum, offset+1, 1, num,
00716 tdata.p0(), &status);
00717 nrows_ = max(nrows_,offset+num);
00718 check_errors();
00719 break;
00720 }
00721 default:
00722 throw Message_error ("unsupported data type in write_column_raw_void()");
00723 }
00724 }
00725
00726 template<typename T> void fitshandle::write_image (const arr2<T> &data)
00727 {
00728 assert_image_hdu("fitshandle::write_image()");
00729 planck_assert (axes_.size()==2, "wrong number of dimensions");
00730 planck_assert (axes_[0]==data.size1(), "wrong size of dimension 1");
00731 planck_assert (axes_[1]==data.size2(), "wrong size of dimension 2");
00732
00733 fits_write_img (fptr, FITSUTIL<T>::DTYPE, 1, axes_[0]*axes_[1],
00734 const_cast<T *>(&data[0][0]), &status);
00735 check_errors();
00736 }
00737
00738 template void fitshandle::write_image (const arr2<float> &data);
00739 template void fitshandle::write_image (const arr2<double> &data);
00740 template void fitshandle::write_image (const arr2<int> &data);
00741
00742 template<typename T> void fitshandle::write_subimage
00743 (const arr<T> &data, int64 offset)
00744 {
00745 assert_image_hdu("fitshandle::write_subimage()");
00746 fits_write_img (fptr, FITSUTIL<T>::DTYPE, 1+offset,
00747 data.size(), const_cast<T *>(&data[0]), &status);
00748 check_errors();
00749 }
00750
00751 template void fitshandle::write_subimage (const arr<float> &data, int64 offset);
00752 template void fitshandle::write_subimage
00753 (const arr<double> &data, int64 offset);
00754 template void fitshandle::write_subimage (const arr<int> &data, int64 offset);
00755
00756 template<typename T> void fitshandle::read_image (arr2<T> &data) const
00757 {
00758 assert_image_hdu("fitshandle::read_image()");
00759 planck_assert (axes_.size()==2, "wrong number of dimensions");
00760 data.alloc(axes_[0], axes_[1]);
00761 fits_read_img (fptr, FITSUTIL<T>::DTYPE, 1, axes_[0]*axes_[1], 0,
00762 &data[0][0], 0, &status);
00763 check_errors();
00764 }
00765
00766 template void fitshandle::read_image (arr2<float> &data) const;
00767 template void fitshandle::read_image (arr2<double> &data) const;
00768 template void fitshandle::read_image (arr2<int> &data) const;
00769
00770 template<typename T> void fitshandle::read_image (arr3<T> &data) const
00771 {
00772 assert_image_hdu("fitshandle::read_image()");
00773 planck_assert (axes_.size()==3, "wrong number of dimensions");
00774 data.alloc(axes_[0], axes_[1], axes_[2]);
00775 fits_read_img (fptr, FITSUTIL<T>::DTYPE, 1, axes_[0]*axes_[1]*axes_[2],
00776 0, &data(0,0,0), 0, &status);
00777 check_errors();
00778 }
00779
00780 template void fitshandle::read_image (arr3<float> &data) const;
00781 template void fitshandle::read_image (arr3<double> &data) const;
00782 template void fitshandle::read_image (arr3<int> &data) const;
00783
00784 template<typename T> void fitshandle::read_subimage
00785 (arr2<T> &data, int xl, int yl) const
00786 {
00787 assert_image_hdu("fitshandle::read_subimage()");
00788 planck_assert (axes_.size()==2, "wrong number of dimensions");
00789 for (int m=0; m<data.size1(); ++m)
00790 fits_read_img (fptr, FITSUTIL<T>::DTYPE, (xl+m)*axes_[1]+yl+1,
00791 data.size2(), 0, &data[m][0], 0, &status);
00792 check_errors();
00793 }
00794
00795 template void fitshandle::read_subimage
00796 (arr2<float> &data, int xl, int yl) const;
00797 template void fitshandle::read_subimage
00798 (arr2<double> &data, int xl, int yl) const;
00799 template void fitshandle::read_subimage
00800 (arr2<int> &data, int xl, int yl) const;
00801
00802 template<typename T> void fitshandle::read_subimage
00803 (arr<T> &data, int64 offset) const
00804 {
00805 assert_image_hdu("fitshandle::read_subimage()");
00806 fits_read_img (fptr, FITSUTIL<T>::DTYPE, 1+offset,
00807 data.size(), 0, &data[0], 0, &status);
00808 check_errors();
00809 }
00810
00811 template void fitshandle::read_subimage (arr<float> &data, int64 offset) const;
00812 template void fitshandle::read_subimage (arr<double> &data, int64 offset) const;
00813 template void fitshandle::read_subimage (arr<int> &data, int64 offset) const;
00814
00815 void fitshandle::add_healpix_keys (int datasize)
00816 {
00817 int nside = isqrt(datasize/12);
00818 planck_assert (12*nside*nside==datasize, "Wrong Healpix map size");
00819
00820 update_key ("PIXTYPE",string("HEALPIX"),"HEALPIX pixelisation");
00821 update_key ("ORDERING",string("RING"),
00822 "Pixel ordering scheme, either RING or NESTED");
00823 update_key ("NSIDE",nside,"Resolution parameter for HEALPIX");
00824 update_key ("FIRSTPIX",0,"First pixel # (0 based)");
00825 update_key ("LASTPIX",datasize-1,"Last pixel # (0 based)");
00826 update_key ("INDXSCHM",string("IMPLICIT"),
00827 "Indexing : IMPLICIT or EXPLICIT");
00828 update_key ("GRAIN",0,"Grain of pixel indexing");
00829 }
00830
00831 namespace {
00832
00833 class cfitsio_checker
00834 {
00835 public:
00836 cfitsio_checker()
00837 {
00838 float fitsversion;
00839 planck_assert(fits_get_version(&fitsversion),
00840 "error calling fits_get_version()");
00841 int v_header = nearest<int>(1000.*CFITSIO_VERSION),
00842 v_library = nearest<int>(1000.*fitsversion);
00843 if (v_header!=v_library)
00844 cerr << endl << "WARNING: version mismatch between CFITSIO header (v"
00845 << dataToString(v_header*0.001) << ") and linked library (v"
00846 << dataToString(v_library*0.001) << ")." << endl << endl;
00847 }
00848 };
00849
00850 static cfitsio_checker Cfitsio_Checker;
00851
00852 }
