// Copyright (C) 1999-2018
// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
// For conditions of distribution and use, see copyright notice in "copyright"

#include <pthread.h>

#include "fitsimage.h"
#include "framebase.h"

#include "context.h"
#include "mmap.h"
#include "smmap.h"
#include "mmapincr.h"
#include "alloc.h"
#include "allocgz.h"
#include "channel.h"
#include "share.h"
#include "sshare.h"
#include "socket.h"
#include "socketgz.h"
#include "var.h"
#include "order.h"
#include "iis.h"
#include "hist.h"
#include "compress.h"
#include "analysis.h"
#include "photo.h"
#include "colorscale.h"

FitsImage::FitsImage(Context* cx, Tcl_Interp* pp)
{
  context_ =cx;
  interp_ =pp;

  objectKeyword_ =NULL;
  fileName =NULL;
  rootBaseFileName =NULL;
  fullBaseFileName =NULL;
  iisFileName =NULL;

  fits_ =NULL;
  post_ =NULL;
  hist_ =NULL;
  hpx_ =NULL;
  base_ =NULL;
  basedata_ =NULL;

  manageBlock_ =0;
  block_ =NULL;
  blockdata_ =NULL;

  manageAnalysis_ =0;
  analysis_ =NULL;
  analysisdata_ =NULL;

  image_ =NULL;
  data_ =NULL;

  nextMosaic_ =NULL;
  nextSlice_ =NULL;

  keyLTMV =0;
  keyATMV =0;
  keyDTMV =0;
  keyDATASEC =0;

  jyBeam_ =1;

  imageToData = Translate(-.5, -.5);
  dataToImage = Translate( .5,  .5);

  imageToData3d = Translate3d(-.5, -.5, -.5);
  dataToImage3d = Translate3d( .5,  .5,  .5);

  imageToRef3d = imageToData3d;
  refToImage3d = dataToImage3d;

  manageWCS_ =1;
  ast_ =NULL;
  astInv_ =0;
  wcs_ =NULL;
  wcsCel_ =NULL;
  wcs3D_ =NULL;
  wcsHPX_ =0;
  wcsSize_ =NULL;

  wcsAltHeader_ =NULL;
  wfpc2Header_ =NULL;
  wcs0Header_ =NULL;

  iisMode_ =0;
  iiszt_ =0;

  for (int ii=0; ii<FTY_MAXAXES; ii++)
    address[ii] =1;
}

FitsImage::~FitsImage()
{
  if (objectKeyword_)
    delete [] objectKeyword_;

  if (fileName)
    delete [] fileName;

  if (rootBaseFileName)
    delete [] rootBaseFileName;

  if (fullBaseFileName)
    delete [] fullBaseFileName;

  if (iisFileName)
    delete [] iisFileName;

  if (fits_)
    delete fits_;

  if (post_)
    delete post_;

  if (hist_)
    delete hist_;

  if (hpx_)
    delete hpx_;

  if (basedata_)
    delete basedata_;

  if (manageBlock_) {
    if (block_)
      delete block_;
    if (blockdata_)
      delete blockdata_;
  }
  if (manageAnalysis_) {
    if (analysis_)
      delete analysis_;
    if (analysisdata_)
      delete analysisdata_;
  }

  if (manageWCS_) {
    if (ast_)
      astAnnul(ast_);
    if (wcs_)
      delete [] wcs_;
    if (wcsCel_)
      delete [] wcsCel_;
    if (wcs3D_)
      delete [] wcs3D_;
    if (wcsSize_)
      delete [] wcsSize_;
  }

  if (wcsAltHeader_)
    delete wcsAltHeader_;
  if (wfpc2Header_)
    delete wfpc2Header_;
  if (wcs0Header_)
    delete wcs0Header_;
}

// Fits

FitsImageFitsAlloc::FitsImageFitsAlloc(Context* cx, Tcl_Interp* pp, 
				       const char* ch, const char* fn, 
				       FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsAlloc(ch, FitsFile::RELAX, flush);
  process(fn,id);
}

FitsImageFitsAllocGZ::FitsImageFitsAllocGZ(Context* cx, Tcl_Interp* pp, 
					   const char* ch, const char* fn, 
					   FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsAllocGZ(ch, FitsFile::RELAX, flush);
  process(fn,id);
}

FitsImageFitsChannel::FitsImageFitsChannel(Context* cx, Tcl_Interp* pp, 
					   const char* ch, const char* fn, 
					   FitsFile::FlushMode flush, int id) 
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsChannel(pp, ch, fn, FitsFile::RELAX, flush);
  process(fn,id);
}

FitsImageFitsMMap::FitsImageFitsMMap(Context* cx, Tcl_Interp* pp, 
				     const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsMMap(fn, FitsFile::RELAX);
  process(fn,id);
}

FitsImageFitsSMMap::FitsImageFitsSMMap(Context* cx, Tcl_Interp* pp, 
				       const char* hdr, 
				       const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsSMMap(hdr, fn);
  process(fn,id);
}

FitsImageFitsMMapIncr::FitsImageFitsMMapIncr(Context* cx, Tcl_Interp* pp, 
					     const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsMMapIncr(fn, FitsFile::RELAX);
  process(fn,id);
}

FitsImageFitsShare::FitsImageFitsShare(Context* cx, Tcl_Interp* pp, 
				       Base::ShmType type, 
				       int sid, const char* fn, int id)
  : FitsImage(cx, pp)
{
  switch (type) {
  case Base::SHMID:
    fits_ = new FitsFitsShareID(sid, fn, FitsFile::RELAX);
    break;
  case Base::KEY:
    fits_ = new FitsFitsShareKey(sid, fn, FitsFile::RELAX);
    break;
  }
  process(fn,id);
}

FitsImageFitsSShare::FitsImageFitsSShare(Context* cx, Tcl_Interp* pp, 
					 Base::ShmType type,
					 int hdr, int sid, 
					 const char* fn, int id)
  : FitsImage(cx, pp)
{
  switch (type) {
  case Base::SHMID:
    fits_ = new FitsFitsSShareID(hdr, sid, fn);
    break;
  case Base::KEY:
    fits_ = new FitsFitsSShareKey(hdr, sid, fn);
    break;
  }
  process(fn,id);
}

FitsImageFitsSocket::FitsImageFitsSocket(Context* cx, Tcl_Interp* pp, 
					 int s, const char* fn, 
					 FitsFile::FlushMode flush, int id) 
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsSocket(s, fn, FitsFile::RELAX, flush);
  process(fn,id);
}

FitsImageFitsSocketGZ::FitsImageFitsSocketGZ(Context* cx, Tcl_Interp* pp,
					     int s, const char* fn, 
					     FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsSocketGZ(s, fn, FitsFile::RELAX, flush);
  process(fn,id);
}

FitsImageFitsVar::FitsImageFitsVar(Context* cx, Tcl_Interp* pp, 
				   const char* var, const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsVar(pp, var, fn, FitsFile::RELAX);
  process(fn,id);
}

FitsImageFitsOrder::FitsImageFitsOrder(Context* cx, Tcl_Interp* pp, 
				       FitsImage* fi,
				       FitsHead* hdr, char* data, size_t sz,
				       int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsOrder(fi->fitsFile(), hdr, data, sz);
  process(NULL,id);

  rootBaseFileName = dupstr(fi->rootBaseFileName);
  fullBaseFileName = dupstr(fi->fullBaseFileName);
  iisFileName = dupstr(fi->fullBaseFileName);
}

// Fits Next

FitsImageFitsNextAlloc::FitsImageFitsNextAlloc(Context* cx, Tcl_Interp* pp, 
					       const char* fn,
					       FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextAlloc(prev);
  process(fn,id);
}

FitsImageFitsNextAllocGZ::FitsImageFitsNextAllocGZ(Context* cx, Tcl_Interp* pp, 
						   const char* fn,
						   FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextAllocGZ(prev);
  process(fn,id);
}

FitsImageFitsNextChannel::FitsImageFitsNextChannel(Context* cx, Tcl_Interp* pp, 
						   const char* fn,
						   FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextChannel(prev);
  process(fn,id);
}

FitsImageFitsNextMMap::FitsImageFitsNextMMap(Context* cx, Tcl_Interp* pp, 
					     const char* fn,
					     FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextMMap(prev);
  process(fn,id);
}

FitsImageFitsNextSMMap::FitsImageFitsNextSMMap(Context* cx, Tcl_Interp* pp, 
					       const char* fn,
					       FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextSMMap(prev);
  process(fn,id);
}

FitsImageFitsNextMMapIncr::FitsImageFitsNextMMapIncr(Context* cx, Tcl_Interp* pp, 
						     const char* fn,
						     FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextMMapIncr(prev);
  process(fn,id);
}

FitsImageFitsNextShare::FitsImageFitsNextShare(Context* cx, Tcl_Interp* pp, 
					       const char* fn,
					       FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextShare(prev);
  process(fn,id);
}

FitsImageFitsNextSShare::FitsImageFitsNextSShare(Context* cx, Tcl_Interp* pp, 
						 const char* fn,
						 FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextSShare(prev);
  process(fn,id);
}

FitsImageFitsNextSocket::FitsImageFitsNextSocket(Context* cx, Tcl_Interp* pp, 
						 const char* fn,
						 FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextSocket(prev);
  process(fn,id);
}

FitsImageFitsNextSocketGZ::FitsImageFitsNextSocketGZ(Context* cx, Tcl_Interp* pp, 
						     const char* fn, 
						     FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextSocketGZ(prev);
  process(fn,id);
}

FitsImageFitsNextVar::FitsImageFitsNextVar(Context* cx, Tcl_Interp* pp, 
					   const char* fn,
					   FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsFitsNextVar(prev);
  process(fn,id);
}

FitsImageFitsNextHist::FitsImageFitsNextHist(Context* cx, Tcl_Interp* pp, 
					     FitsImage* fi,
					     FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsHistNext(prev);
  process(NULL,id);

  fits_->setpFilter(fi->getHistFilter());
  fits_->setpBinX(fi->getHistX());
  fits_->setpBinY(fi->getHistY());
  fits_->setpBinZ(fi->getHistZ());

  rootBaseFileName = dupstr(fi->rootBaseFileName);
  fullBaseFileName = dupstr(fi->fullBaseFileName);
  iisFileName = dupstr(fi->fullBaseFileName);
}

FitsImageFitsNextPost::FitsImageFitsNextPost(Context* cx, Tcl_Interp* pp, 
					     FitsImage* fi,
					     FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsPostNext(prev);
  process(NULL,id);

  rootBaseFileName = dupstr(fi->rootBaseFileName);
  fullBaseFileName = dupstr(fi->fullBaseFileName);
  iisFileName = dupstr(fi->fullBaseFileName);
}

FitsImageFitsNextOrder::FitsImageFitsNextOrder(Context* cx, Tcl_Interp* pp, 
					       FitsImage* fi,
					       FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsOrderNext(prev);
  process(NULL,id);

  rootBaseFileName = dupstr(fi->rootBaseFileName);
  fullBaseFileName = dupstr(fi->fullBaseFileName);
  iisFileName = dupstr(fi->fullBaseFileName);
}

// Array

FitsImageArrAlloc::FitsImageArrAlloc(Context* cx, Tcl_Interp* pp, 
				     const char* ch, const char* fn,
				     FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrAlloc(ch, flush);
  process(fn,id);
}

FitsImageArrAllocGZ::FitsImageArrAllocGZ(Context* cx, Tcl_Interp* pp, 
					 const char* ch, const char* fn,
					 FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrAllocGZ(ch, flush);
  process(fn,id);
}

FitsImageArrChannel::FitsImageArrChannel(Context* cx, Tcl_Interp* pp, 
					 const char* ch, const char* fn,
					 FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrChannel(pp, ch, fn, flush);
  process(fn,id);
}

FitsImageArrMMap::FitsImageArrMMap(Context* cx, Tcl_Interp* pp, 
				   const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrMMap(fn);
  process(fn,id);
}

FitsImageArrMMapIncr::FitsImageArrMMapIncr(Context* cx, Tcl_Interp* pp, 
					   const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrMMapIncr(fn);
  process(fn,id);
}

FitsImageArrShare::FitsImageArrShare(Context* cx, Tcl_Interp* pp, 
				     Base::ShmType type, 
				     int sid, const char* fn, int id) 
  : FitsImage(cx, pp)
{
  switch (type) {
  case Base::SHMID:
    fits_ = new FitsArrShareID(sid, fn);
    break;
  case Base::KEY:
    fits_ = new FitsArrShareKey(sid, fn);
    break;
  }
  process(fn,id);
}

FitsImageArrSocket::FitsImageArrSocket(Context* cx, Tcl_Interp* pp, 
				       int s, const char* fn,
				       FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrSocket(s, fn, flush);
  process(fn,id);
}

FitsImageArrSocketGZ::FitsImageArrSocketGZ(Context* cx, Tcl_Interp* pp, 
					   int s, const char* fn,
					   FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrSocketGZ(s, fn, flush);
  process(fn,id);
}

FitsImageArrVar::FitsImageArrVar(Context* cx, Tcl_Interp* pp, 
				 const char* var, const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsArrVar(pp, var, fn);
  process(fn,id);
}

// ENVI

FitsImageENVISMMap::FitsImageENVISMMap(Context* cx, Tcl_Interp* pp,
				       const char* hdr,
				       const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsENVISMMap(hdr,fn);
  process(fn,id);
}

// NRRD

FitsImageNRRDAlloc::FitsImageNRRDAlloc(Context* cx, Tcl_Interp* pp,
				     const char* ch, const char* fn,
				     FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsNRRDAlloc(ch, flush);
  process(fn,id);
}

FitsImageNRRDChannel::FitsImageNRRDChannel(Context* cx, Tcl_Interp* pp,
					 const char* ch, const char* fn,
					 FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsNRRDChannel(pp, ch, fn, flush);
  process(fn,id);
}

FitsImageNRRDMMap::FitsImageNRRDMMap(Context* cx, Tcl_Interp* pp,
				     const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsNRRDMMap(fn);
  process(fn,id);
}

FitsImageNRRDShare::FitsImageNRRDShare(Context* cx, Tcl_Interp* pp,
				       Base::ShmType type, 
				       int sid, const char* fn, int id) 
  : FitsImage(cx, pp)
{
  switch (type) {
  case Base::SHMID:
    fits_ = new FitsNRRDShareID(sid, fn);
    break;
  case Base::KEY:
    fits_ = new FitsNRRDShareKey(sid, fn);
    break;
  }
  process(fn,id);
}

FitsImageNRRDSocket::FitsImageNRRDSocket(Context* cx, Tcl_Interp* pp,
					 int s, const char* fn,
					 FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsNRRDSocket(s, fn, flush);
  process(fn,id);
}

FitsImageNRRDVar::FitsImageNRRDVar(Context* cx, Tcl_Interp* pp,
				   const char* var, const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsNRRDVar(pp, var, fn);
  process(fn,id);
}

// Photo

FitsImagePhoto::FitsImagePhoto(Context* cx, Tcl_Interp* pp,
			       const char* ph, const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsPhoto(pp, ph);
  process(fn,id);
}

FitsImagePhotoCube::FitsImagePhotoCube(Context* cx, Tcl_Interp* pp,
				       const char* ph, const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsPhotoCube(pp, ph);
  process(fn,id);
}

FitsImagePhotoCubeNext::FitsImagePhotoCubeNext(Context* cx, Tcl_Interp* pp,
					       const char* fn,
					       FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsPhotoCubeNext(prev);
  process(fn,id);
}

// Mosaic

FitsImageMosaicAlloc::FitsImageMosaicAlloc(Context* cx, Tcl_Interp* pp,
					   const char* ch, 
					   const char* fn,
					   FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicAlloc(ch, flush);
  process(fn,id);
}

FitsImageMosaicAllocGZ::FitsImageMosaicAllocGZ(Context* cx, Tcl_Interp* pp,
					       const char* ch, 
					       const char* fn,
					       FitsFile::FlushMode flush, 
					       int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicAllocGZ(ch, flush);
  process(fn,id);
}

FitsImageMosaicChannel::FitsImageMosaicChannel(Context* cx, Tcl_Interp* pp,
					       const char* ch, 
					       const char* fn,
					       FitsFile::FlushMode flush,
					       int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicChannel(pp, ch, flush);
  process(fn,id);
}

FitsImageMosaicMMap::FitsImageMosaicMMap(Context* cx, Tcl_Interp* pp,
					 const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicMMap(fn);
  process(fn,id);
}

FitsImageMosaicMMapIncr::FitsImageMosaicMMapIncr(Context* cx, Tcl_Interp* pp,
						 const char* fn, 
						 int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicMMapIncr(fn);
  process(fn,id);
}

FitsImageMosaicShare::FitsImageMosaicShare(Context* cx, Tcl_Interp* pp,
					   Base::ShmType type,
					   int sid, const char* fn, int id)
  : FitsImage(cx, pp)
{
  switch (type) {
  case Base::SHMID:
    fits_ = new FitsMosaicShareID(sid);
    break;
  case Base::KEY:
    fits_ = new FitsMosaicShareKey(sid);
    break;
  }
  process(fn,id);
}

FitsImageMosaicSocket::FitsImageMosaicSocket(Context* cx, Tcl_Interp* pp,
					     int s, const char* fn,
					     FitsFile::FlushMode flush, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicSocket(s, flush);
  process(fn,id);
}

FitsImageMosaicSocketGZ::FitsImageMosaicSocketGZ(Context* cx, Tcl_Interp* pp,
						 int s, const char* fn,
						 FitsFile::FlushMode flush, 
						 int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicSocketGZ(s, flush);
  process(fn,id);
}

FitsImageMosaicVar::FitsImageMosaicVar(Context* cx, Tcl_Interp* pp,
				       const char* var, const char* fn, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicVar(pp, var, fn);
  process(fn,id);
}

// Mosaic Next

FitsImageMosaicNextAlloc::FitsImageMosaicNextAlloc(Context* cx, Tcl_Interp* pp,
						   const char* fn,
						   FitsFile* prev,
						   FitsFile::FlushMode flush,
						   int id)

  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextAlloc(prev, flush);
  process(fn,id);
}

FitsImageMosaicNextAllocGZ::FitsImageMosaicNextAllocGZ(Context* cx, Tcl_Interp* pp,
						       const char* fn,
						       FitsFile* prev,
						      FitsFile::FlushMode flush,
						       int id)

  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextAllocGZ(prev, flush);
  process(fn,id);
}

FitsImageMosaicNextChannel::FitsImageMosaicNextChannel(Context* cx, 
						       Tcl_Interp* pp,
						       const char* fn, 
						       FitsFile* prev,
						      FitsFile::FlushMode flush,
						       int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextChannel(prev, flush);
  process(fn,id);
}

FitsImageMosaicNextMMap::FitsImageMosaicNextMMap(Context* cx, Tcl_Interp* pp,
						 const char* fn, 
						 FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextMMap(prev);
  process(fn,id);
}

FitsImageMosaicNextMMapIncr::FitsImageMosaicNextMMapIncr(Context* cx, Tcl_Interp* pp,
							 const char* fn, 
							 FitsFile* prev,
							 int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextMMapIncr(prev);
  process(fn,id);
}

FitsImageMosaicNextShare::FitsImageMosaicNextShare(Context* cx, Tcl_Interp* pp, 
						   const char* fn,
						   FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextShare(prev);
  process(fn,id);
}

FitsImageMosaicNextSocket::FitsImageMosaicNextSocket(Context* cx, 
						     Tcl_Interp* pp,
						     const char* fn, 
						     FitsFile* prev,
						     FitsFile::FlushMode flush,
						     int id)

  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextSocket(prev, flush);
  process(fn,id);
}

FitsImageMosaicNextSocketGZ::FitsImageMosaicNextSocketGZ(Context* cx, 
							 Tcl_Interp* pp,
							 const char* fn, 
							 FitsFile* prev,
						      FitsFile::FlushMode flush,
							 int id)

  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextSocketGZ(prev, flush);
  process(fn,id);
}

FitsImageMosaicNextVar::FitsImageMosaicNextVar(Context* cx, Tcl_Interp* pp,
					       const char* fn, 
					       FitsFile* prev, int id)
  : FitsImage(cx, pp)
{
  fits_ = new FitsMosaicNextVar(prev);
  process(fn,id);
}

// IIS

FitsImageIIS::FitsImageIIS(Context* cx, Tcl_Interp* pp, int w, int h)
  : FitsImage(cx, pp)
{
  fits_ = new FitsIIS(w, h);
  process("",1);
  iisMode_ = 1;
}

void FitsImageIIS::iisErase()
{
  ((FitsIIS*)fits_)->erase();
}

char* FitsImageIIS::iisGet(int xx, int yy, int dx, int dy)
{
  return ((FitsIIS*)fits_)->get(xx, yy, dx, dy);
}

void FitsImageIIS::iisSet(const char* src, int xx, int yy, int dx, int dy)
{
  ((FitsIIS*)fits_)->set(src, xx, yy, dx, dy);
}

void FitsImageIIS::iisWCS(const Matrix& mm, const Vector& z, int zt)
{
  Matrix& mx = (Matrix&)mm;
  double sx = mx[0][0];
  double sy = mx[1][1];
  double tx = mx[2][0];
  double ty = mx[2][1];

  imageToPhysical = 
    Translate(0,-height()/2.) *
    FlipY() * 
    Translate(0,height()/2.) *
    Translate(-1,0) * 
    Scale(sx,sy) * 
    Translate(tx,ty);
  physicalToImage = imageToPhysical.invert();
  iisz_ = z;
  iiszt_ = zt;
}

// FitsImage

void FitsImage::wfpc2WCS(istream& str)
{
  FitsHead* hh = parseWCS(str);

  // Process OBJECT keyword
  if (objectKeyword_)
    delete [] objectKeyword_;
  objectKeyword_ = dupstr(hh->getString("OBJECT"));

  // Process WCS keywords
  if (wfpc2Header_)
    delete wfpc2Header_;

  wfpc2Header_ = hh;
  initWCS(wfpc2Header_);
}

void FitsImage::appendWCS(istream& str)
{
  FitsHead* hh = parseWCS(str);

  // process OBJECT keyword
  char* obj = dupstr(hh->getString("OBJECT"));
  if (obj) {
    if (objectKeyword_)
      delete [] objectKeyword_;
    objectKeyword_ = obj;
  }

  // Process WCS keywords
  FitsHead* hd = image_->head();

  // append wcs keywords to the end of the header
  int ll = hd->headbytes()+hh->headbytes();
  char* cards = new char[ll];

  // copy default wcs
  memcpy(cards, hd->cards(), hd->headbytes());

  // find first END and zero out
  for (int i=0; i<hd->headbytes(); i+=80)
    if (!strncmp(cards+i,"END",3)) {
      memcpy(cards+i, "   ",3);
      break;
    }

  // copy appended wcs
  memcpy(cards+hd->headbytes(), hh->cards(), hh->headbytes());

  delete hh;
  if (wcsAltHeader_)
    delete wcsAltHeader_;

  wcsAltHeader_ = new FitsHead(cards,ll,FitsHead::ALLOC);
  initWCS(wcsAltHeader_);
}

char* FitsImage::display(FitsHead* hd)
{
  int size = hd->ncard() * (FTY_CARDLEN+1);
  char* lbuf = new char[size+1];

  char* lptr = lbuf;
  char* cptr = hd->cards();
  for (int i=0; i<hd->ncard(); i++,cptr+=FTY_CARDLEN) {
    memcpy(lptr, cptr, FTY_CARDLEN);
    lptr+=FTY_CARDLEN;
    *(lptr++) = '\n';
  }

  lbuf[size] = '\0';
  return lbuf;
}

char* FitsImage::displayHeader()
{
  Vector blockFactor = context_->blockFactor();
  if (blockFactor[0] != 1 && blockFactor[1] != 1)
    return display(image_->head());

  if (DebugBin || DebugCompress)
    return display(image_->head());
  else
    return display(fits_->head());
}

char* FitsImage::displayPrimary()
{
  return display(fits_->primary());
}

char* FitsImage::displayWCS()
{
  if (wcsAltHeader_)
    return display(wcsAltHeader_);
  else if (wfpc2Header_)
    return display(wfpc2Header_);
  else if (wcs0Header_)
    return display(wcs0Header_);
  else
    return display(image_->head());
}

FitsBound* FitsImage::getDataParams(FrScale::SecMode which)
{
  switch (which) {
  case FrScale::IMGSEC:
    return &iparams;
  case FrScale::DATASEC:
    return &dparams;
  case FrScale::CROPSEC:
    return &cparams;
  }

  // just to satisfy the compiler
  return &iparams;
}

const char* FitsImage::getValue(const Vector& v)
{
  if (!isIIS())
    return data_->getValue(v);
  else {
    double val = data_->getValueDouble(v);

    ostringstream str;
    if (val == 0)
      str << ends;
    else if (val == IISMIN)
      str << '<' << iisz_[0] << ends;
    else if (val == IISMAX)
      str << '>' << iisz_[1] << ends;
    else if (val > IISMAX)
      str << ends;
    else
      // W_LINEAR =1
      if (iiszt_ == 1)
	str << ((val-IISMIN) * (iisz_[1]-iisz_[0]))/(IISMAX-IISMIN) + iisz_[0] 
	    << ends;
      else
	str << val << ends;

    memcpy(buf,str.str().c_str(), str.str().length());
    return buf;
  }
}

void FitsImage::iisSetFileName(const char* fn)
{
  if (iisFileName)
    delete [] iisFileName;
  iisFileName = dupstr(fn);
}

void FitsImage::initWCS(FitsHead* hd)
{
  if (manageWCS_) {
    if (ast_)
      astAnnul(ast_);
    ast_ =NULL;
    if (wcs_)
      delete [] wcs_;
    wcs_ =NULL;
    if (wcs_)
      delete [] wcs_;
    wcs_ =NULL;
    if (wcsCel_)
      delete [] wcsCel_;
    wcsCel_ =NULL;
    if (wcs3D_)
      delete [] wcs3D_;
    wcs3D_ =NULL;
    wcsHPX_ = 0;
    if (wcsSize_)
      delete [] wcsSize_;
    wcsSize_ =NULL;
  }

  // shareWCS?
  manageWCS_ =1;
  if (context_->shareWCS()) {
    FitsImage* ptr = context_->fits;
    while (ptr) {
      if (ptr == this)
	break;

      FitsImage* sptr = ptr->nextSlice();
      while (sptr) {
	if (sptr == this) {
	  ast_ = ptr->ast_;
	  astInv_ = ptr->astInv_;
	  wcs_ = ptr->wcs_;
	  wcsCel_ = ptr->wcsCel_;
	  wcs3D_ = ptr->wcs3D_;
	  wcsHPX_ = ptr->wcsHPX_;
	  wcsSize_ = ptr->wcsSize_;

	  initWCSPhysical();
	  manageWCS_ =0;
	  return;
	}
	sptr = sptr->nextSlice();
      }
      ptr = ptr->nextMosaic();
    }
  }

  int hasWCSAST = hd->find("BEGAST_A") ? 1 : 0;

  astInit(hd);
  wcsInit(hasWCSAST);
  wcsCelInit(hasWCSAST);
  wcs3DInit(hasWCSAST);
  wcsHPXInit();
  wcsSizeInit();
  
  initWCSPhysical();

  if (DebugWCS && ast_)
    astShow(ast_);
}

void FitsImage::resetWCS()
{
  // Process OBJECT keyword
  if (objectKeyword_)
    delete [] objectKeyword_;
  objectKeyword_ = dupstr(image_->getString("OBJECT"));

  if (wcsAltHeader_)
    delete wcsAltHeader_;
  wcsAltHeader_ =NULL;

  if (wcs0Header_)
    delete wcs0Header_;
  wcs0Header_ =NULL;

  if (wfpc2Header_)
    initWCS(wfpc2Header_);
  else
    initWCS(image_->head());
}

void FitsImage::initWCS0(const Vector& pix)
{
  if (!ast_)
    return;

  FitsHead* hd = new FitsHead(naxis(0), naxis(1), 1, -32);

  // CTYPE
  hd->appendString("CTYPE1", "RA---TAN", NULL);
  hd->appendString("CTYPE2", "DEC--TAN", NULL);

  // CRPIX
  Vector cc = mapFromRef(pix, Coord::IMAGE, Coord::FK5);
  hd->appendReal("CRPIX1", cc[1], 8, NULL);
  hd->appendReal("CRPIX2", cc[0], 8, NULL);

  // CRVAL
  hd->appendReal("CRVAL1", 0, 8, NULL);
  hd->appendReal("CRVAL2", 0, 8, NULL);

  // CD
  double ss = getWCSSize(Coord::WCS);
  double ang = getWCSRotation(Coord::WCS,Coord::FK5);
  Matrix flip;
  switch (getWCSOrientation(Coord::WCS,Coord::FK5)) {
  case Coord::NORMAL:
  case Coord::YY:
    flip = FlipX();
    break;
  case Coord::XX:
  case Coord::XY:
    break;
  };
  Matrix mx = flip*Rotate(ang)*Scale(ss);
  hd->appendReal("CD1_1", mx[0][0], 8, NULL);
  hd->appendReal("CD1_2", mx[0][1], 8, NULL);
  hd->appendReal("CD2_1", mx[1][0], 8, NULL);
  hd->appendReal("CD2_2", mx[1][1], 8, NULL);

  // EPOCH, EQUINOX
  hd->appendReal("EPOCH", 2000, 8, NULL);
  hd->appendReal("EQUINOX", 2000, 8, NULL);
  
  // RADESYS
  hd->appendString("RADESYS", "FK5", NULL);

  if (wcs0Header_)
    delete wcs0Header_;
  wcs0Header_ = hd;
  initWCS(wcs0Header_);
}

void FitsImage::initWCSPhysical()
{
  // now see if we have a 'physical' in WCSP, if so, set LTMV keywords
  keyLTMV =0;

  char* wcsname = image_->getString("WCSNAMEP");
  if (wcsname && *wcsname && !strncmp(wcsname, "PHYSICAL", 8)) {
    if (image_->find("CD1_1P")  || image_->find("CD1_2P")  ||
	image_->find("CD2_1P")  || image_->find("CD2_2P")  ||
	image_->find("CRPIX1P") || image_->find("CRPIX2P") ||
	image_->find("CRVAL1P") || image_->find("CRVAL2P")) {
      keyLTMV = 1;

      double cd11 = image_->getReal("CD1_1P", 1);
      double cd12 = image_->getReal("CD1_2P", 0);
      double cd21 = image_->getReal("CD2_1P", 0);
      double cd22 = image_->getReal("CD2_2P", 1);

      double crpix1 = image_->getReal("CRPIX1P", 0);
      double crpix2 = image_->getReal("CRPIX2P", 0);
      double crval1 = image_->getReal("CRVAL1P", 0);
      double crval2 = image_->getReal("CRVAL2P", 0);

      double ltm11 = cd11 != 0 ? 1/cd11 : 0;
      double ltm12 = cd12 != 0 ? 1/cd12 : 0;
      double ltm21 = cd21 != 0 ? 1/cd21 : 0;
      double ltm22 = cd22 != 0 ? 1/cd22 : 0;

      double ltv1 = crpix1 - crval1*ltm11 - crval2*ltm21;
      double ltv2 = crpix2 - crval1*ltm12 - crval2*ltm22;

      physicalToImage = Matrix(ltm11, ltm12, ltm21, ltm22, ltv1, ltv2);
      imageToPhysical = physicalToImage.invert();
    }
  }
}

void FitsImage::load()
{
  if (post_)
    base_ = post_;
  else if (hpx_)
    base_ = hpx_;
  else if (hist_)
    base_ = hist_;
  else
    base_ = fits_;

  if (basedata_)
    delete basedata_;

  switch (base_->head()->bitpix()) {
  case 8:
    basedata_ = new FitsDatam<unsigned char>(base_, interp_);
    break;
  case 16:
    basedata_ = new FitsDatam<short>(base_, interp_);
    break;
  case -16:
    basedata_ = new FitsDatam<unsigned short>(base_, interp_);
    break;
  case 32:
    basedata_ = new FitsDatam<int>(base_, interp_);
    break;
  case 64:
    basedata_ = new FitsDatam<long long>(base_, interp_);
    break;
  case -32:
    basedata_ = new FitsDatam<float>(base_, interp_);
    break;
  case -64:
    basedata_ = new FitsDatam<double>(base_, interp_);
    break;
  }

  // block
  block_ = base_;
  blockdata_ = basedata_;

  // analysis
  analysis_ = block_;
  analysisdata_ = blockdata_;

  // image
  image_ = analysis_;
  data_ = analysisdata_;
}

void FitsImage::match(const char* xxname1, const char* yyname1,
		      Coord::CoordSystem sys1, Coord::SkyFrame sky1,
		      const char* xxname2, const char* yyname2,
		      Coord::CoordSystem sys2, Coord::SkyFrame sky2,
		      double rad, Coord::CoordSystem sys, 
		      Coord::DistFormat dist,
		      const char* rrname)
{
  if (!hasWCS(sys1) || !hasWCS(sys2))
    return;

  // only good for skyframe

  astClearStatus; // just to make sure
  astBegin; // start memory management

  // get lists
  Tcl_Obj* listxx1 =
    Tcl_GetVar2Ex(interp_, xxname1, NULL, TCL_LEAVE_ERR_MSG);
  Tcl_Obj* listyy1 =
    Tcl_GetVar2Ex(interp_, yyname1, NULL, TCL_LEAVE_ERR_MSG);
  Tcl_Obj* listxx2 =
    Tcl_GetVar2Ex(interp_, xxname2, NULL, TCL_LEAVE_ERR_MSG);
  Tcl_Obj* listyy2 =
    Tcl_GetVar2Ex(interp_, yyname2, NULL, TCL_LEAVE_ERR_MSG);

  // get objects
  int nxx1;
  Tcl_Obj **objxx1;
  Tcl_ListObjGetElements(interp_, listxx1, &nxx1, &objxx1);
  int nyy1;
  Tcl_Obj **objyy1;
  Tcl_ListObjGetElements(interp_, listyy1, &nyy1, &objyy1);
  int nxx2;
  Tcl_Obj **objxx2;
  Tcl_ListObjGetElements(interp_, listxx2, &nxx2, &objxx2);
  int nyy2;
  Tcl_Obj **objyy2;
  Tcl_ListObjGetElements(interp_, listyy2, &nyy2, &objyy2);

  // sanity check
  if (nxx1 != nyy1 || nxx2 != nyy2)
    return;

  if (!hasWCSCel(sys1) || !hasWCSCel(sys2))
    return;

  // get doubles
  double* ixx1 = new double[nxx1];
  for (int ii=0 ; ii<nxx1 ; ii++)
    Tcl_GetDoubleFromObj(interp_, objxx1[ii], ixx1+ii);
  double* iyy1 = new double[nyy1];
  for (int ii=0 ; ii<nyy1 ; ii++)
    Tcl_GetDoubleFromObj(interp_, objyy1[ii], iyy1+ii);
  double* ixx2 = new double[nxx2];
  for (int ii=0 ; ii<nxx2 ; ii++)
    Tcl_GetDoubleFromObj(interp_, objxx2[ii], ixx2+ii);
  double* iyy2 = new double[nyy2];
  for (int ii=0 ; ii<nyy2 ; ii++)
    Tcl_GetDoubleFromObj(interp_, objyy2[ii], iyy2+ii);

  Vector* in1 = new Vector[nxx1];
  for (int ii=0; ii<nxx1; ii++)
    in1[ii] = degToRad(Vector(ixx1[ii],iyy1[ii]));

  Vector* ptr2 = new Vector[nxx2];
  for (int ii=0; ii<nxx2; ii++)
    ptr2[ii] = degToRad(Vector(ixx2[ii],iyy2[ii]));

  double rr = rad;
  switch (dist) {
  case Coord::DEGREE:
    break;
  case Coord::ARCMIN:
    rr /= 60.;
    break;
  case Coord::ARCSEC:
    rr /= 60.*60.;
    break;
  }
  rr = zeroTWOPI(degToRad(rr));

  Vector* ptr1 =NULL;
  if (sky1 != sky2) {
    setWCSSkyFrame(sys1, sky1);
    AstFrameSet* wcs1 = (AstFrameSet*)astCopy(ast_);

    setWCSSkyFrame(sys2, sky2);
    AstFrameSet* wcs2 = (AstFrameSet*)astCopy(ast_);

    AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, "SKY");
    if (cvt != AST__NULL) {
      ptr1 = new Vector[nxx1];
      wcsTran(cvt, nxx1, in1, 1, ptr1);
    }
  }
  else
    ptr1 = in1;

  // now compare
  if (ptr1 && ptr2) {
    setWCSSkyFrame(sys2, sky2);
    Tcl_Obj* objrr = Tcl_NewListObj(0,NULL);
    for(int jj=0; jj<nxx2; jj++) {
      for (int ii=0; ii<nxx1; ii++) {
	double dd = wcsDistance(ptr1[ii], ptr2[jj]);

	if ((dd != AST__BAD) && (dd <= rr)) {
	  Tcl_Obj* obj[2];
	  obj[0] = Tcl_NewIntObj(ii+1);
	  obj[1] = Tcl_NewIntObj(jj+1);
	  Tcl_Obj* list = Tcl_NewListObj(2,obj);
	  Tcl_ListObjAppendElement(interp_, objrr, list);
	}
      }
    }
    Tcl_SetVar2Ex(interp_, rrname, NULL, objrr, TCL_LEAVE_ERR_MSG);
  }

  // clean up
  astEnd; // now, clean up memory

  // clean up
  if (ixx1)
    delete [] ixx1;
  if (iyy1)
    delete [] iyy1;
  if (ixx2)
    delete [] ixx2;
  if (iyy2)
    delete [] iyy2;

  if (!ptr1 && ptr1!=in1)
    delete [] ptr1;
  if (in1)
    delete [] in1;
  if (ptr2)
    delete [] ptr2;
}

Matrix& FitsImage::matrixToData(Coord::InternalSystem sys)
{
  switch (sys) {
  case Coord::WINDOW:
    return windowToData;
  case Coord::CANVAS:
    return canvasToData;
  case Coord::WIDGET:
    return widgetToData;
  case Coord::PANNER:
    return pannerToData;
  case Coord::MAGNIFIER:
    return magnifierToData;
  case Coord::PS:
    return psToData;
  case Coord::USER:
    return userToData;
  case Coord::REF:
    return refToData;
  }

  // just to satisfy the compiler
  return refToData;
}

Matrix& FitsImage::matrixFromData(Coord::InternalSystem sys)
{
  switch (sys) {
  case Coord::WINDOW:
    return dataToWindow;
  case Coord::CANVAS:
    return dataToCanvas;
  case Coord::WIDGET:
    return dataToWidget;
  case Coord::PANNER:
    return dataToPanner;
  case Coord::MAGNIFIER:
    return dataToMagnifier;
  case Coord::PS:
    return dataToPS;
  case Coord::USER:
    return dataToUser;
  case Coord::REF:
    return dataToRef;
  }

  // just to satisfy the compiler
  return dataToRef;
}

Matrix3d& FitsImage::matrixFromData3d(Coord::InternalSystem sys)
{
  switch (sys) {
  case Coord::WINDOW:
    return dataToWindow3d;
  case Coord::CANVAS:
    return dataToCanvas3d;
  case Coord::WIDGET:
    return dataToWidget3d;
  case Coord::PANNER:
    return dataToPanner3d;
  case Coord::MAGNIFIER:
    return dataToMagnifier3d;
  case Coord::PS:
    return dataToPS3d;
  case Coord::USER:
    return dataToUser3d;
  case Coord::REF:
    return dataToRef3d;
  }

  // just to satisfy the compiler
  return dataToRef3d;
}

Matrix3d& FitsImage::matrixToData3d(Coord::InternalSystem sys)
{
  switch (sys) {
  case Coord::WINDOW:
    return windowToData3d;
  case Coord::CANVAS:
    return canvasToData3d;
  case Coord::WIDGET:
    return widgetToData3d;
  case Coord::PANNER:
    return pannerToData3d;
  case Coord::MAGNIFIER:
    return magnifierToData3d;
  case Coord::PS:
    return psToData3d;
  case Coord::USER:
    return userToData3d;
  case Coord::REF:
    return refToData3d;
  }

  // just to satisfy the compiler
  return refToData3d;
}

FitsHead* FitsImage::parseWCS(istream& str)
{
  FitsHead* hd = image_->head();
  FitsHead* rr = new FitsHead(hd->naxis(0), hd->naxis(1),
			      hd->naxis(2), hd->bitpix());

  // this works for both terminated (\n) and non-terminated (fits) headers
  while (!str.eof()) {
    char buf[256];
    str.get(buf,80,'\n');
    if (strlen(buf) > 0) {
      // check for blank lines
      if (*buf == ' ')
	break;

      string x(buf);
      istringstream sstr(x);

      char keyword[32];
      sstr >> keyword;

      if (strchr(buf,'=')) {
	char dummy;
	sstr >> dummy;
      }

      if (strchr(buf,'\'')) {
	char val[64];
	char* ss = strchr(buf,'\'')+1;
	char* ee = strrchr(buf,'\'');
	int ll = ee-ss;
	if (ll<0 || ll>63)
	  ll =0;
	strncpy(val,ss,ll);
	val[ll] = '\0';
	rr->appendString(keyword, val, "");
      }
      else {
	double val;
	sstr >> val;
	rr->appendReal(keyword, val, 15, "");
      }

      if (strlen(buf) <= 80) {
	// eat the \n
	char b;
	str.get(b);
      }
    }
    else
      break;
  }

  return rr;
}

void FitsImage::process(const char* fn, int id)
{
  if (!fits_->isValid()) {
    reset();
    return;
  }

  if (fits_->isImage()) {
    switch (fits_->pEncoding()) {
    case FitsFile::RAW:
    case FitsFile::BSQ:
      break;
    case FitsFile::GZIP:
      initNRRD();
      if (!post_ || !post_->isValid()) {
	reset();
	return;
      }
      break;
    case FitsFile::BIL:
    case FitsFile::BIP:
      initENVI();
      if (!post_ || !post_->isValid()) {
	reset();
	return;
      }
      break;
    default:
      reset();
      return;
    }

    load();
  }
  else if (fits_->isBinTable()) {
    // Compress
    if (fits_->find("ZIMAGE")) {
      initCompress();
      if (!post_ || !post_->isValid()) {
	reset();
	return;
      }
      load();
    }
    // HEALPIX
    else if ((fits_->find("PIXTYPE") && (!strncmp(fits_->getString("PIXTYPE"),"HEALPIX",4))) || fits_->find("NSIDE")) {
      initHPX();
      if (!hpx_ || !hpx_->isValid()) {
	reset();
	return;
      }
      load();
    }
    else {
      // Bintable
      initBin();
      if (!hist_ || !hist_->isValid()) {
	reset();
	return;
      }
    }
  }
  else if (fits_->isAsciiTable()) {
    // HEALPIX
    if (fits_->find("NSIDE")) {
      initHPX();
      if (!hpx_ || !hpx_->isValid()) {
	reset();
	return;
      }
      load();
    }
  }

  // set slice address
  for (int ii=1; ii<id; ii++) {
    for (int jj=2; jj<FTY_MAXAXES; jj++) {
      if (address[jj]<naxis(jj)) {
	address[jj]++;
	break;
      }
      else
	address[jj]=1;
    }
  }

  // load() can call reset()
  if (fits_)
    setFileName(fn);
}

void FitsImage::processKeywordsPhysical()
{
  // Physical to Image (LTM/LTV keywords) (with no wcsname already located)

  if (!keyLTMV) {
    if (image_->find("LTM1_1") ||
	image_->find("LTM1_2") ||
	image_->find("LTM2_1") ||
	image_->find("LTM2_2") ||
	image_->find("LTV1") ||
	image_->find("LTV2"))
      keyLTMV = 1;

    double ltm11 = image_->getReal("LTM1_1", 1);
    double ltm12 = image_->getReal("LTM1_2", 0);
    double ltm21 = image_->getReal("LTM2_1", 0);
    double ltm22 = image_->getReal("LTM2_2", 1);

    double ltv1 = image_->getReal("LTV1", 0);
    double ltv2 = image_->getReal("LTV2", 0);

    physicalToImage = Matrix(ltm11, ltm12, ltm21, ltm22, ltv1, ltv2);
    imageToPhysical = physicalToImage.invert();
  }

  // CDD to Detector (DTM/DTV keywords)

  keyDTMV =0;
  if (image_->find("DTM1_1") ||
      image_->find("DTM1_2") ||
      image_->find("DTM2_1") ||
      image_->find("DTM2_2") ||
      image_->find("DTV1") ||
      image_->find("DTV2"))
    keyDTMV = 1;

  double dtm11 = image_->getReal("DTM1_1", 1);
  double dtm12 = image_->getReal("DTM1_2", 0);
  double dtm21 = image_->getReal("DTM2_1", 0);
  double dtm22 = image_->getReal("DTM2_2", 1);

  double dtv1 = image_->getReal("DTV1", 0);
  double dtv2 = image_->getReal("DTV2", 0);

  physicalToDetector = Matrix(dtm11, dtm12, dtm21, dtm22, dtv1, dtv2);
  detectorToPhysical = physicalToDetector.invert();

  // Physical to Amplifier (ATM/ATV keywords)

  keyATMV =0;
  if (image_->find("ATM1_1") ||
      image_->find("ATM1_2") ||
      image_->find("ATM2_1") ||
      image_->find("ATM2_2") ||
      image_->find("ATV1") ||
      image_->find("ATV2"))
    keyATMV = 1;

  double atm11 = image_->getReal("ATM1_1", 1);
  double atm12 = image_->getReal("ATM1_2", 0);
  double atm21 = image_->getReal("ATM2_1", 0);
  double atm22 = image_->getReal("ATM2_2", 1);

  double atv1 = image_->getReal("ATV1", 0);
  double atv2 = image_->getReal("ATV2", 0);

  physicalToAmplifier = Matrix(atm11, atm12, atm21, atm22, atv1, atv2);
  amplifierToPhysical = physicalToAmplifier.invert();

  if (DebugMosaic) {
    cerr << endl;
    cerr << "ATM/V: " << physicalToAmplifier << endl;
    cerr << "ATM/V-1: " << amplifierToPhysical << endl;
    cerr << "DTM/V: " << physicalToDetector << endl;
    cerr << "DTM/V-1: " << detectorToPhysical << endl;
    cerr << "LTM/V: " << physicalToImage << endl;
    cerr << "LTM/V-1: " << imageToPhysical << endl;
  }

  /*
  // Radio data?
  char* bunit = image_->getString("BUNIT");
  double cdelt1 = fabs(image_->getReal("CDELT1",0));
  double cdelt2 = fabs(image_->getReal("CDELT2",0));
  double bmaj = image_->getReal("BMAJ",0);
  double bmin = image_->getReal("BMIN",0);

  // ok we have a radio map
#define GFACTOR (2.0*sqrt(2.0*log(2.0)))
  if (!strncmp(bunit,"JY/BEAM",7) && cdelt1 && cdelt2 && bmaj && bmin) {
    // convert from deg to arcsec?
    cdelt1 *= 3600;
    cdelt2 *= 3600;
    bmaj *= 3600;
    bmin *= 3600;
    jyBeam_ = (2*M_PI*bmaj*bmin)/(GFACTOR*GFACTOR*cdelt1*cdelt2);
  }

  if (bunit)
    delete [] bunit;
  if (cunit1)
    delete [] cunit1;
  if (cunit2)
    delete [] cunit2;
  */
}

void FitsImage::processKeywordsParams()
{
  // iparams is a BBOX in DATA coords 0-n
  iparams.set(0, 0, width(), height());

  {
    char* datstr = image_->getString("DATASEC");
    // default
    Vector v1(1,1);
    Vector v2(size());

    keyDATASEC =0;
    if (datstr && *datstr && parseSection(datstr,&v1,&v2)) {
      // additional check
      if (v1[0]<1 || v1[1]<1 || 
	  v1[1]>width() || v2[1]>height() ||
	  v1[0]>v2[0] || v1[1]>v2[1]) {
	// default
	v1 = Vector(1,1);
	v2 = Vector(size());
	keyDATASEC = 0;
      } else
	keyDATASEC = 1;
    }

    // dparams is a BBOX in DATA coords 0-n
    datasec = BBox(v1,v2);
    v1 -= Vector(1,1);
    dparams.set(v1[0],v1[1],v2[0],v2[1]);
  }

  // DEBUG
  if (DebugCrop) {
    cerr << "iparams " << iparams << endl;
    cerr << "dparams " << dparams << endl;
  }
}

void FitsImage::processKeywordsFitsSection()
{
  Vector ll(dparams.xmin,dparams.ymin);
  Vector ur(dparams.xmax,dparams.ymax);

  if (fits_->pcoord() && fits_->pxvalid() && fits_->pyvalid()) {
    ll[0] = fits_->pxmin();
    ur[0] = fits_->pxmax();
    ll[1] = fits_->pymin();
    ur[1] = fits_->pymax();

    ll = ll*physicalToImage*Translate(-1,-1);
    ur = ur*physicalToImage;

    context_->setSecMode(FrScale::CROPSEC);
  }
  if (!fits_->pcoord() && fits_->pxvalid()) {
    ll[0] = fits_->pxmin()-1;
    ur[0] = fits_->pxmax();
    context_->setSecMode(FrScale::CROPSEC);
  }
  if (!fits_->pcoord() && fits_->pyvalid()) {
    ll[1] = fits_->pymin()-1;
    ur[1] = fits_->pymax();
    context_->setSecMode(FrScale::CROPSEC);
  }

  // params is a BBOX in DATA coords 0-n
  setCropParams(ll,ur,0);

  // DEBUG
  if (DebugCrop)
    cerr << "cparams " << cparams << endl;

  if (fits_->pzvalid()) {
    int zmin = fits_->pzmin()-1;
    int zmax = fits_->pzmax();
    context_->setSecMode(FrScale::CROPSEC);
    context_->setCrop3dParams(zmin,zmax);
  }
}

int FitsImage::processKeywordsIRAF(FitsImage* fits)
{
  // DETSEC
  Coord::Orientation orientation = Coord::NORMAL;

  char* detstr =  image_->getString("DETSEC");
  Vector dv1,dv2;
  if (!(detstr && *detstr && parseSection(detstr,&dv1,&dv2)))
    return 0;

  BBox detsec = BBox(dv1,dv2);

  int xx = (dv1[0] < dv2[0]);
  int yy = (dv1[1] < dv2[1]);

  if (xx && yy)
    orientation = Coord::NORMAL;
  else if (!xx && yy)
    orientation = Coord::XX;
  else if (!xx && !yy)
    orientation = Coord::XY;
  else if (xx && !yy)
    orientation = Coord::YY;

  // DETSIZE
  char* sizestr = image_->getString("DETSIZE");
  Vector sv1(1,1);
  Vector sv2(10000,10000);
  if (sizestr && *sizestr) {
    if (!(parseSection(sizestr,&sv1,&sv2)))
      return 0;
  }

  BBox detsize = BBox(sv1,sv2);
  
  // CCDSUM
  Vector ccdsum(1,1);
  char* ccdstr = image_->getString("CCDSUM");
  if (ccdstr && *ccdstr) {
    double Ns, Np, Ns1, Np1;
    string x(ccdstr);
    istringstream str(x);

    str >> Ns >> Np >> Ns1 >> Np1;
    ccdsum = Vector(1/Ns, 1/Np);
  }
  
  // origin
  Vector origin = detsec.ll * Scale(ccdsum) * Translate(-datasec.ll);

  // matrix
  // if the segment is flipped, we can have a discontinuity at
  // the edges, due to round off errors, so we 'nudge' it

  Matrix flip;
  switch (orientation) {
  case Coord::NORMAL:
    break;
  case Coord::XX:
    flip = FlipX();
      break;
  case Coord::YY:
    flip = FlipY();
      break;
  case Coord::XY:
    flip = FlipXY();
      break;
  }

  // internal flip
  Matrix mflip;
  switch (context_->IRAFOrientation(orientation)) {
  case Coord::NORMAL:
    break;
  case Coord::XX:
    mflip = FlipX();
    break;
  case Coord::YY:
    mflip = FlipY();
    break;
  case Coord::XY:
    mflip = FlipXY();
    break;
  }

  Vector center = datasec.center() * imageToData;
  Vector mcenter = detsize.center() * imageToData * Scale(ccdsum);

  wcsToRef_ =
    Translate(-center) *
    flip *
    Translate(center) *
    Translate(origin) * 
    Translate(-mcenter) *
    mflip *
    Translate(mcenter);

  // we do this to shift the origin to the middle of the image to match
  // the wcs case. Needed by imageBBox()
  // first? reset wcsToRef
  if (fits == this) {
    irafToRef = wcsToRef_.invert();
    wcsToRef_ = Matrix();
  }
  else
    wcsToRef_ *= fits->irafToRef;

  if (DebugMosaic) {
    cerr << "ProcessKeywordsIRAF" << endl
	 << " datasec: " << datasec << endl
	 << " ccdsum : " << ccdsum << endl
	 << " detsize: " << detsize << endl
	 << " detsec : " << detsec << endl
	 << " matrix : " << wcsToRef_ << endl;
  }

  return 1;
}

void FitsImage::replaceWCS(istream& str)
{
  FitsHead* hh = parseWCS(str);

  // Process OBJECT keyword
  if (objectKeyword_)
    delete [] objectKeyword_;
  objectKeyword_ = dupstr(hh->getString("OBJECT"));

  // Process WCS keywords
  if (wcsAltHeader_)
    delete wcsAltHeader_;

  wcsAltHeader_ = hh;
  initWCS(wcsAltHeader_);
}

void FitsImage::reset()
{
  if (fits_)
    delete fits_;
  fits_ =NULL;

  if (post_)
    delete post_;
  post_ =NULL;
  if (hpx_)
    delete hpx_;
  hpx_ =NULL;
  if (hist_)
    delete hist_;
  hist_ =NULL;

  base_ =NULL;
  if (basedata_)
    delete basedata_;
  basedata_ =NULL;

  if (manageBlock_) {
    if (block_)
      delete block_;
    if (blockdata_)
      delete blockdata_;
  }
  manageBlock_ =0;
  block_ =NULL;
  blockdata_ =NULL;

  if (manageAnalysis_) {
    if (analysis_)
      delete analysis_;
    if (analysisdata_)
      delete analysisdata_;
  }
  manageAnalysis_ =0;
  analysis_ =NULL;
  analysisdata_ =NULL;

  image_ =NULL;
  data_ =NULL;
}

char* FitsImage::root(const char* fn)
{
  if (fn) {
    const char* ptr = fn;           // init the ptr
    while(*ptr++);                  // walk it forward to end of string
    ptr--;                          // backup one
    while(*ptr != '/' && ptr != fn) // walk it backward til last / or beginning
      ptr--;
    if (*ptr == '/')                // step it over the last '/'
      ptr++;
    return dupstr(ptr);             // got it!
  }
  else
    return NULL;
}

void FitsImage::setCropParams(int datasec)
{
  if (!datasec)
    cparams = iparams;
  else
    cparams = dparams;
}

void FitsImage::setCropParams(const Vector& ss, const Vector& tt, int datasec)
{
  // Coord are in DATA
  Vector ll = ss;
  Vector ur = tt;

  int xmin = ll[0];
  int xmax = ur[0];
  int ymin = ll[1];
  int ymax = ur[1];

  if (xmin>xmax) {
    xmin = ur[0];
    xmax = ll[0];
  }
  if (ymin>ymax) {
    ymin = ur[1];
    ymax = ll[1];
  }

  setCropParams(xmin,ymin,xmax,ymax,datasec);
}

void FitsImage::setCropParams(int x0, int y0, int x1, int y1, int datasec)
{
  FitsBound* params;
  if (!datasec)
    params = &iparams;
  else
    params = &dparams;

  // Coords are in DATA
  if (x0<params->xmin)
    x0=params->xmin;
  if (x0>params->xmax)
    x0=params->xmax;
  if (x1<params->xmin)
    x1=params->xmin;
  if (x1>params->xmax)
    x1=params->xmax;

  if (y0<params->ymin)
    y0=params->ymin;
  if (y0>params->ymax)
    y0=params->ymax;
  if (y1<params->ymin)
    y1=params->ymin;
  if (y1>params->ymax)
    y1=params->ymax;

  cparams.set(x0,y0,x1,y1);
}

void FitsImage::setFileName(const char* fn)
{
  if (fileName)
    delete [] fileName;
  fileName = NULL;

  if (rootBaseFileName)
    delete [] rootBaseFileName;
  rootBaseFileName = NULL;

  if (fullBaseFileName)
    delete [] fullBaseFileName;
  fullBaseFileName = NULL;

  if (iisFileName)
    delete [] iisFileName;
  iisFileName = NULL;

  // no filename to set
  if (!fn)
    return;

  // strip any '[]'
  char* ffn = strip(fn);

  FitsFile* ptr = post_ ? post_ : fits_;
  if (!ptr)
    return;

  const char* ext = ptr->extname();
  if (ext) {
    {
      ostringstream str;
      str << ffn << '[' << ext << ']' << ends;
      fullBaseFileName = dupstr(str.str().c_str());
    }
    {
      char* m = root(ffn);
      ostringstream str;
      str << m << '[' << ext << ']' << ends;
      rootBaseFileName = dupstr(str.str().c_str());
      delete [] m;
    }
  }
  else if (ptr->ext()) {
    {
      ostringstream str;
      str << ffn << '[' << ptr->ext() << ']' << ends;
      fullBaseFileName = dupstr(str.str().c_str());
    }
    {
      char* m = root(ffn);
      ostringstream str;
      str << m << '[' << ptr->ext() << ']' << ends;
      rootBaseFileName = dupstr(str.str().c_str());
      delete [] m;
    }
  }
  else {
    fullBaseFileName = dupstr(ffn);
    rootBaseFileName = root(ffn);
  }
  
  // by default, iisFileName is fullBaseFileName
  if (fullBaseFileName)
    iisFileName = dupstr(fullBaseFileName);

  delete [] ffn;
}

void FitsImage::setObjectKeyword(const char* obj)
{
  if (objectKeyword_)
    delete [] objectKeyword_;
  objectKeyword_ = dupstr(obj);
}

char* FitsImage::strip(const char* fn)
{
  if (fn) {
    char* r = dupstr(fn);           // dup the string
    char* ptr = r;                  // init the ptr
    while(*ptr != '[' && *ptr)      // walk it forward til '[' or end
      ptr++;
    *ptr = '\0';                    // zero out rest

    return r;                       // got it!
  }
  else
    return NULL;
}

int FitsImage::nhdu()
{
  int dd =1;
  for (int ii=2; ii<FTY_MAXAXES; ii++)
    if (naxis(ii))
      dd *= naxis(ii);
  return dd;
}

void FitsImage::updateClip(FrScale* fr)
{
  data_->updateClip(fr,getDataParams(fr->secMode()));
}

void* clipproc(void* tt)
{
  t_clip_arg* targ = (t_clip_arg*)tt;
  FitsData* data = targ->data;
  FrScale* fr = targ->fr;
  FitsBound* bb = targ->bb;
  data->updateClip(fr,bb);

  return NULL;
}

void FitsImage::updateClip(FrScale* fr, pthread_t* thread, t_clip_arg* targ)
{
  targ->data = data_;
  targ->fr = fr;
  targ->bb = getDataParams(fr->secMode());
  int result = pthread_create(thread, NULL, clipproc, targ);
  if (result)
    internalError("Unable to Create Thread");
}

const char* FitsImage::getFileName(Base::FileNameType type)
{
  switch (type) {
  case Base::ROOTBASE:
    return rootBaseFileName;
  case Base::FULLBASE:
    return fullBaseFileName;
  case Base::ROOT:
  case Base::FULL:
    // clear the buffer
    if (fileName)
      delete [] fileName;
    fileName =NULL;

    // generate filename
    // if FITS bin table cube, be sure to check the first slice
    if (context_->fits->isHist())
      return updateFileNameBin(type);
    else
      return updateFileNameImage(type);
  }

  // just to satisfy the compiler
  return rootBaseFileName;
}

const char* FitsImage::updateFileNameImage(Base::FileNameType type)
{
  // 2d/3d section
  char* sec =NULL;
  switch (context_->secMode()) {
  case FrScale::IMGSEC:
  case FrScale::DATASEC:
    {
      Vector blockFactor = context_->blockFactor();

      if (blockFactor[0] != 1) {
	ostringstream str;
	str << "*," << blockFactor[0] << ends;
	sec = dupstr(str.str().c_str());
      }
    }    
    break;
  case FrScale::CROPSEC:
    {
      FitsBound* params =getDataParams(FrScale::CROPSEC);

      // params is a BBOX in DATA coords 0-n
      //   xlate to 1-n
      Vector ll= Vector(params->xmin,params->ymin) * Translate(1,1);
      Vector ur(params->xmax,params->ymax);
      Vector blockFactor = context_->blockFactor();

      if (blockFactor[0] != 1) {
	ostringstream str;
	str << ll[0] << ':' << ur[0] << ',' 
	    << ll[1] << ':' << ur[1] << ','
	    << blockFactor[0] << ends;
	sec = dupstr(str.str().c_str());
      }
      else {
	ostringstream str;
	str << ll[0] << ':' << ur[0] << ',' 
	    << ll[1] << ':' << ur[1] << ends;
	sec = dupstr(str.str().c_str());
      }
    }    
    break;
  }

  // address
  char* add =NULL;
  {
    int doAdd =0;
    ostringstream str;

    int jj;
    for (jj=FTY_MAXAXES-1; jj>=2; jj--) {
      if (address[jj]!=1)
	break;
    }	
    jj++;
    for (int ii=2; ii<jj; ii++) {
      doAdd =1;
      if (ii==2)
	str << "plane=";
      else
	str << ':';
      str << address[ii];
    }
    if (doAdd) {
      str << ends;
      add = dupstr(str.str().c_str());
    }
  }

  switch (type) {
  case Base::ROOTBASE:
  case Base::FULLBASE:
    // better not get here
    return NULL;

  case Base::ROOT:
    if (rootBaseFileName) {
      ostringstream str;
      if (add && sec)
	str << rootBaseFileName << '[' << add << ']' << '[' << sec << ']';
      else if (add && !sec)
	str << rootBaseFileName << '[' << add << ']';
      else if (!add && sec)
	str << rootBaseFileName << '[' << sec << ']';
      else
	str << rootBaseFileName;
      str << ends;
      fileName = dupstr(str.str().c_str());
    }
    break;

  case Base::FULL:
    if (fullBaseFileName) {
      ostringstream str;
      if (add && sec)
	str << fullBaseFileName << '[' << add << ']' << '[' << sec << ']';
      else if (add && !sec)
	str << fullBaseFileName << '[' << add << ']';
      else if (!add && sec)
	str << fullBaseFileName << '[' << sec << ']';
      else
	str << fullBaseFileName;
      str << ends;
      fileName = dupstr(str.str().c_str());
    }
    break;
  }

  if (sec)
    delete [] sec;
  if (add)
    delete [] add;

  return fileName;
}

const char* FitsImage::updateFileNameBin(Base::FileNameType type)
{
  char* filter = (char*)fits_->pFilter();
  int doFilter = (filter && *filter);
    
  // x,y filter
  char* sec = NULL;
  {
    switch (context_->secMode()) {
    case FrScale::IMGSEC:
    case FrScale::DATASEC:
      break;
    case FrScale::CROPSEC:
      {
	ostringstream str;
	FitsBound* params =getDataParams(FrScale::CROPSEC);

	// dataToPhysical not set yet
	Vector ll = Vector(params->xmin,params->ymin) *
	  dataToImage * imageToPhysical;
	Vector ur = Vector(params->xmax,params->ymax) * 
	  dataToImage * imageToPhysical;

	str << fits_->pBinX() << ">=" << ll[0] << ',' 
	    << fits_->pBinX() << "<=" << ur[0] << ','
	    << fits_->pBinY() << ">=" << ll[1] << ',' 
	    << fits_->pBinY() << "<=" << ur[1] << ends;
	sec = dupstr(str.str().c_str());
      }    
      break;
    }
  }

  // z filter
  char* slice =NULL;
  {
    char* zcol = (char*)fits_->pBinZ();
    int bd = context_->binDepth();
    if (bd>1 && zcol) {
      // only the first slice will have this
      FitsImage* first = context_->fits;
      if (first) {
	Vector zlim = first->fits_->getColMinMax(zcol);
	double zlen = zlim[1]-zlim[0];
	double zdelta = zlen/bd;

	double zptr = zlim[0] + (address[2]-1)*zdelta;
	ostringstream str;
	str << zcol << ">=" << zptr << '&' << zcol << '<' << zptr+zdelta <<ends;
	slice = dupstr(str.str().c_str());
      }
    }
  }

  switch (type) {
  case Base::ROOTBASE:
  case Base::FULLBASE:
    // better not get here
    return NULL;

  case Base::ROOT:
    if (rootBaseFileName) {
      ostringstream str;
      str << rootBaseFileName;
      if (doFilter && sec && slice)
	str << '[' << filter << ',' << sec << ',' << slice << ']';
      else if (doFilter && sec && !slice)
	str << '[' << filter << ',' << sec << ']';
      else if (doFilter && !sec && slice)
	str << '[' << filter << ',' << slice << ']';
      else if (doFilter && !sec && !slice)
	str << '[' << filter << ']';
      else if (!doFilter && sec && slice)
	str << '[' << sec << ',' << slice << ']';
      else if (!doFilter && sec && !slice)
	str << '[' << sec << ']';
      else if (!doFilter && !sec && slice)
	str << '[' << slice << ']';
      str << ends;
      fileName = dupstr(str.str().c_str());
    }
  case Base::FULL:
    if (fullBaseFileName) {
      ostringstream str;
      str << fullBaseFileName;
      if (doFilter && sec && slice)
	str << '[' << filter << ',' << sec << ',' << slice << ']';
      else if (doFilter && sec && !slice)
	str << '[' << filter << ',' << sec << ']';
      else if (doFilter && !sec && slice)
	str << '[' << filter << ',' << slice << ']';
      else if (doFilter && !sec && !slice)
	str << '[' << filter << ']';
      else if (!doFilter && sec && slice)
	str << '[' << sec << ',' << slice << ']';
      else if (!doFilter && sec && !slice)
	str << '[' << sec << ']';
      else if (!doFilter && !sec && slice)
	str << '[' << slice << ']';
      str << ends;
      fileName = dupstr(str.str().c_str());
    }
    break;
  }

  return fileName;
}

void FitsImage::updateMatrices(Matrix& rgbToRef,
			       Matrix& refToUser,
			       Matrix& userToWidget,
			       Matrix& widgetToCanvas,
			       Matrix& canvasToWindow)
{
  dataToRef = wcsToRef_ * rgbToRef;
  refToData = dataToRef.invert();

  dataToUser = dataToRef * refToUser;
  userToData = dataToUser.invert();

  dataToWidget = dataToUser * userToWidget;
  widgetToData = dataToWidget.invert();

  dataToCanvas = dataToWidget * widgetToCanvas;
  canvasToData = dataToCanvas.invert();

  dataToWindow = dataToCanvas * canvasToWindow;
  windowToData = dataToWindow.invert();

  refToCanvas = refToUser * userToWidget * widgetToCanvas;
  canvasToRef = refToCanvas.invert();

  imageToRef = imageToData * dataToRef;
  refToImage = imageToRef.invert();

  imageToWidget = imageToRef * refToUser * userToWidget;
  widgetToImage = imageToWidget.invert();

  physicalToRef = physicalToImage * imageToData * dataToRef;
  refToPhysical = physicalToRef.invert();

  amplifierToRef = amplifierToPhysical * physicalToRef;
  refToAmplifier = amplifierToRef.invert();

  detectorToRef = detectorToPhysical * physicalToRef;
  refToDetector = detectorToRef.invert();
}

void FitsImage::updateMatrices(Matrix3d& refToUser3d, 
			       Matrix3d& userToWidget3d,
			       Matrix3d& widgetToCanvas3d, 
			       Matrix3d& canvasToWindow3d)
{
  dataToUser3d = dataToRef3d * refToUser3d;
  userToData3d = dataToUser3d.invert();

  dataToWidget3d = dataToUser3d * userToWidget3d;
  widgetToData3d = dataToWidget3d.invert();

  dataToCanvas3d = dataToWidget3d * widgetToCanvas3d;
  canvasToData3d = dataToCanvas3d.invert();

  dataToWindow3d = dataToCanvas3d * canvasToWindow3d;
  windowToData3d = dataToWindow3d.invert();
}

void FitsImage::updatePannerMatrices(Matrix& refToPanner)
{
  dataToPanner = dataToRef * refToPanner;
  pannerToData = dataToPanner.invert();
}

void FitsImage::updatePannerMatrices(Matrix3d& refToPanner3d)
{
  dataToPanner3d = dataToRef3d * refToPanner3d;
  pannerToData3d = dataToPanner3d.invert();
}

void FitsImage::updateMagnifierMatrices(Matrix& refToMagnifier)
{
  dataToMagnifier = dataToRef * refToMagnifier;
  magnifierToData = dataToMagnifier.invert();
}

void FitsImage::updateMagnifierMatrices(Matrix3d& refToMagnifier3d)
{
  dataToMagnifier3d = dataToRef3d * refToMagnifier3d;
  magnifierToData3d = dataToMagnifier3d.invert();
}

void FitsImage::updatePS(Matrix ps)
{
  dataToPS = dataToRef * ps;
  psToData = dataToPS.invert();
}

void FitsImage::updatePS(Matrix3d ps)
{
  dataToPS3d = dataToRef3d * ps;
  psToData3d = dataToPS3d.invert();
}

// WCS

double FitsImage::getWCSSize(Coord::CoordSystem sys)
{
  if (!wcsSize_ || sys<Coord::WCS)
    return 0;
  else
    return wcsSize_[sys-Coord::WCS];
}

double FitsImage::calcWCSSize(Coord::CoordSystem sys)
{
  if (!hasWCS(sys))
    return 0;

  astClearStatus; // just to make sure
  setWCSSkyFrame(sys, Coord::FK5);
  
  // only check lon axis as it will not change near poles
  Vector in[3];
  Vector out[3];
  in[0] = center();
  in[1] = center()+Vector(0,1);
  wcsTran(2, in, 1, out);

  double dd = wcsDistance(out[0],out[1]);
  return hasWCSCel(sys) ? radToDeg(dd) : dd;
}

Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
						Coord::SkyFrame sky)
{
  if (!hasWCS(sys))
    return Coord::NORMAL;
  
  setWCSSkyFrame(sys, sky);

  Vector in[3];
  Vector out[3];
  in[0] = center();
  in[1] = center()+Vector(0,1);
  in[2] = center()+Vector(1,0);
  wcsTran(3, in, 1, out);
  double ang = wcsAngle(out[1],out[0],out[2]);

  if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX))) {
    if ((hasWCSCel(sys) && ang>0) || (!hasWCSCel(sys) && ang<0))
      return Coord::XX;
    else
      return Coord::NORMAL;
  }
  else
    return Coord::NORMAL;
}

double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
{
  if (!hasWCS(sys))
    return 0;
  
  setWCSSkyFrame(sys, sky);

  if (!wcsHPX_) {
    Vector in[3];
    Vector out[3];
    in[0] = center();
    in[1] = center()+Vector(0,1);
    in[2] = center()+Vector(1,0);
    wcsTran(3, in, 1, out);
    double ang = wcsAxAngle(out[0], out[1]);
    double ang3 = wcsAngle(out[1],out[0],out[2]);

    if ((!(isnan(ang )||isinf(ang) ||(ang ==-DBL_MAX)||(ang ==DBL_MAX))) &&
	(!(isnan(ang3)||isinf(ang3)||(ang3==-DBL_MAX)||(ang3==DBL_MAX)))) {
      if ((hasWCSCel(sys) && ang3>0) || (!hasWCSCel(sys) && ang3<0))
	return -ang;
      else
	return ang;
    }
    else
      return 0;
  }
  else { // special case for HPX
    Vector cc = center();
    Vector wcc = wcsTran(cc, 1);
    Vector wnorth = wcc + Vector(0,.001);
    Vector north = wcsTran(wnorth,0);

    int current = astGetI(ast_,"Current");
    int base = astGetI(ast_,"Base");
    astSetI(ast_,"Current",base);
    double ang = wcsAxAngle(cc,north);
    astSetI(ast_,"Current",current);
    if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX)))
      return ang;
    else
      return 0;
  }
}

const char* FitsImage::getWCSName(Coord::CoordSystem sys) 
{
  if (fits_->find("WCSNAME"))
    return fits_->getString("WCSNAME");
  else
    return NULL;
}

Vector FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
			  Coord::SkyFrame sky)
{
  astClearStatus; // just to make sure

  if (!hasWCS(sys))
    return Vector();
  
  setWCSSkyFrame(sys, sky);

  Vector out = wcsTran(in, 1);
  if (astOK && checkWCS(out))
    return hasWCSCel(sys) ? zero360(radToDeg(out)) : out;
  else
    return Vector();
}

char* FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
			 Coord::SkyFrame sky, Coord::SkyFormat format,
			 char* lbuf)
{
  astClearStatus; // just to make sure
  lbuf[0] = '\0';

  if (!hasWCS(sys))
    return lbuf;

  setWCSSkyFrame(sys, sky);
  
  ostringstream str;
  Vector out = wcsTran(in, 1);
  if (astOK && checkWCS(out)) {
    if (hasWCSCel(sys)) {
      ostringstream hms;
      hms << "hms." << context_->parent_->precHMS_;
      ostringstream dms;
      dms << "+dms." << context_->parent_->precDMS_;

      switch (format) {
      case Coord::DEGREES:
	out = zero360(radToDeg(out));
	str << setprecision(context_->parent_->precDeg_)
	    << out[0] << ' ' << out[1] << ' '
	    << (hasWCSCel(sys) ? coord.skyFrameStr(sky) : "") << ends;
	
	break;

      case Coord::SEXAGESIMAL:
	out = zeroTWOPI(out);
	switch (sky) {
	case Coord::FK4:
	case Coord::FK4_NO_E:
	case Coord::FK5:
	case Coord::ICRS:
	  setWCSFormat(1,hms.str().c_str());
	  setWCSFormat(2,dms.str().c_str());
	  break;
	case Coord::GALACTIC:
	case Coord::SUPERGALACTIC:
	case Coord::ECLIPTIC:
	case Coord::HELIOECLIPTIC:
	  setWCSFormat(1,dms.str().c_str());
	  setWCSFormat(2,dms.str().c_str());
	  break;
	}
	str << astFormat(ast_,1,out[0]) << ' '
	    << astFormat(ast_,2,out[1]) << ' '
	    << (hasWCSCel(sys) ? coord.skyFrameStr(sky) : "") << ends;
	break;
      }
    }
    else
      str << setprecision(context_->parent_->precLinear_)
	  << out[0] << ' ' << out[1] << ends;

    strncpy(lbuf, str.str().c_str(), str.str().length());
  }
  
  return lbuf;
}

Vector3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
			    Coord::SkyFrame sky)
{
  astClearStatus; // just to make sure

  if (!hasWCS(sys))
    return Vector();

  setWCSSkyFrame(sys, sky);

  Vector3d out = wcsTran(in, 1);
  if (astOK && checkWCS(out))
    return hasWCSCel(sys) ? zero360(radToDeg(out)) : out;
  else
    return Vector3d();
}

char* FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
			 Coord::SkyFrame sky, Coord::SkyFormat format,
			 char* lbuf)
{
  astClearStatus; // just to make sure
  lbuf[0] = '\0';

  if (!hasWCS(sys))
    return lbuf;

  setWCSSkyFrame(sys, sky);
  
  ostringstream str;
  Vector3d out = wcsTran(in, 1);
  if (astOK && checkWCS(out)) {
    if (hasWCSCel(sys)) {
      ostringstream hms;
      hms << "hms." << context_->parent_->precHMS_;
      ostringstream dms;
      dms << "+dms." << context_->parent_->precDMS_;

      switch (format) {
      case Coord::DEGREES:
	out = zero360(radToDeg(out));
	str << setprecision(context_->parent_->precDeg_)
	    << out[0] << ' ' << out[1] << ' ' << out[2]
	    << ' ' << (hasWCSCel(sys) ? coord.skyFrameStr(sky) : "") << ends;
	break;

      case Coord::SEXAGESIMAL:
	out = zeroTWOPI(out);
	switch (sky) {
	case Coord::FK4:
	case Coord::FK4_NO_E:
	case Coord::FK5:
	case Coord::ICRS:
	  setWCSFormat(1,hms.str().c_str());
	  setWCSFormat(2,dms.str().c_str());
	  break;
	case Coord::GALACTIC:
	case Coord::SUPERGALACTIC:
	case Coord::ECLIPTIC:
	case Coord::HELIOECLIPTIC:
	  setWCSFormat(1,dms.str().c_str());
	  setWCSFormat(2,dms.str().c_str());
	  break;
	}
	str << astFormat(ast_,1,out[0]) << ' '
	    << astFormat(ast_,2,out[1]) << ' ' << out[2] << ' '
	    << (hasWCSCel(sys) ? coord.skyFrameStr(sky) : "") << ends;
	break;
      }
    }
    else
      str << setprecision(context_->parent_->precLinear_)
	  << out[0] << ' ' << out[1] << ' ' << out[2] <<ends;

    strncpy(lbuf, str.str().c_str(), str.str().length());
  }
  
  return lbuf;
}

Vector FitsImage::wcs2pix(const Vector& vv, Coord::CoordSystem sys,
			  Coord::SkyFrame sky)
{
  astClearStatus; // just to make sure

  if (hasWCS(sys) && astInv_) {
    setWCSSkyFrame(sys, sky);

    Vector in = hasWCSCel(sys) ? degToRad(vv) : vv;
    Vector out = wcsTran(in, 0);
    if (astOK && checkWCS(out))
      return out;
  }
  
  maperr =1;
  return Vector();
}

Vector3d FitsImage::wcs2pix(const Vector3d& vv, Coord::CoordSystem sys,
			    Coord::SkyFrame sky)
{
  astClearStatus; // just to make sure

  if (hasWCS(sys) && astInv_) {
    setWCSSkyFrame(sys, sky);

    Vector3d in = hasWCSCel(sys) ? degToRad(vv) : vv;
    Vector3d out = wcsTran(in, 0);
    if (astOK && checkWCS(out))
      return out;
  }

  return Vector3d();
}

double FitsImage::getWCSDist(const Vector& vv1, const Vector& vv2,
			     Coord::CoordSystem sys)
{
  if (!hasWCS(sys))
    return 0;

  astClearStatus; // just to make sure
  setWCSSkyFrame(sys, Coord::FK5);

  return hasWCSCel(sys) ?
    radToDeg(wcsDistance(degToRad(vv1), degToRad(vv2))) :
    wcsDistance(vv1, vv2);
}

int FitsImage::hasWCS(Coord::CoordSystem sys)
{
  if (!wcs_ || sys<Coord::WCS)
    return 0;
  else
    return wcs_[sys-Coord::WCS];
}

int FitsImage::hasWCSCel(Coord::CoordSystem sys)
{
  if (!wcsCel_ || sys<Coord::WCS)
    return 0;
  else
    return wcsCel_[sys-Coord::WCS];
}

int FitsImage::hasWCSLinear(Coord::CoordSystem sys)
{
  if (!wcs_ || !wcsCel_ || sys<Coord::WCS)
    return 0;
  else
    return wcs_[sys-Coord::WCS] && !wcsCel_[sys-Coord::WCS];
}

int FitsImage::hasWCS3D(Coord::CoordSystem sys)
{
  if (!wcs3D_ || sys<Coord::WCS)
    return 0;
  else
    return wcs3D_[sys-Coord::WCS];
}

void FitsImage::astInit(FitsHead* hd)
{
  if (ast_)
    astAnnul(ast_);
  ast_ =NULL;
  
  // just in case
  if (!hd)
    return;

  ast_ = fits2ast(hd);
  if (!ast_)
    return;

  int naxes = astGetI(ast_,"Naxes");
  switch (naxes) {
  case 1:
    break;
  case 2:
    if (astIsASkyFrame(astGetFrame(ast_,AST__CURRENT)) &&
	astGetI(ast_,"LatAxis") == 1) {
      int orr[] = {2,1};
      astPermAxes(ast_,orr);
    }
    break;
  case 3:
    break;
  case 4:
    break;
  }
}

void FitsImage::wcsInit(int hasWCSAST)
{
  // init wcs_ array
  if (wcs_)
    delete [] wcs_;
  wcs_ =NULL;

  wcs_ = new int[MULTWCS];
  for (int ii=0; ii<MULTWCS; ii++)
    wcs_[ii] =0;
  
  if (!ast_)
    return;

  // since we have ast_
  wcs_[0] =1; 

  // do we have a AST wcs?
  if (hasWCSAST)
    return;

  // fill out wcs_ array
  astClearStatus;
  astBegin;

  int nn = astGetI(ast_, "Nframe");
  for (int ii=0; ii<nn; ii++) {
    const char* id = astGetC(astGetFrame(ast_,ii+1), "Ident");
    if (id && *id) {
      int jj = (*id == ' ') ? 0 : *id-'@';
      wcs_[jj] = 1;
    }
  }

  astEnd;
}

void FitsImage::wcsCelInit(int hasWCSAST)
{
  // init wcsCel_ array
  if (wcsCel_)
    delete [] wcsCel_;
  wcsCel_ =NULL;

  wcsCel_ = new int[MULTWCS];
  for (int ii=0; ii<MULTWCS; ii++)
    wcsCel_[ii] =0;

  if (!ast_)
    return;

  astClearStatus;
  astBegin;

  int nn = astGetI(ast_, "Nframe");
  // do we have a AST wcs?
  if (hasWCSAST) {
    for (int ii=0; ii<nn; ii++) {
      AstFrame* ff = (AstFrame*)astGetFrame(ast_,ii+1);
      AstFrameSet* fs =
	(AstFrameSet*)astFindFrame(ff, astSkyFrame(" MaxAxes=10")," ");
      if (fs) {
	wcsCel_[0] = 1;
	const char* str = astGetC(ff, "System");
	if (!strncmp(str,"Unknown",7))
	  wcsCel_[0] = 0;
      }
    }
  }
  else {
    for (int ii=0; ii<nn; ii++) {
      AstFrame* ff = (AstFrame*)astGetFrame(ast_,ii+1);
      const char* id = astGetC(ff, "Ident");
      if (id && *id) {
	int jj = (*id == ' ') ? 0 : *id-'@';
	AstFrameSet* fs =
	  (AstFrameSet*)astFindFrame(ff, astSkyFrame(" MaxAxes=10")," ");
	if (fs) {
	  wcsCel_[jj] = 1;
	  const char* str = astGetC(ff, "System");
	  if (!strncmp(str,"Unknown",7))
	    wcsCel_[jj] = 0;
	}
      }
    }
  }

  astEnd;
}

void FitsImage::wcs3DInit(int hasWCSAST)
{
  // init wcs3D_ array
  if (wcs3D_)
    delete [] wcs3D_;
  wcs3D_ =NULL;

  wcs3D_ = new int[MULTWCS];
  for (int ii=0; ii<MULTWCS; ii++)
    wcs3D_[ii] =0;
  
  if (!ast_)
    return;

  astClearStatus;
  astBegin;

  int nn = astGetI(ast_,"nframe");
  // do we have a AST wcs?
  if (hasWCSAST) {
    for (int ii=0; ii<nn; ii++) {
      AstFrame* ff = (AstFrame*)astGetFrame(ast_,ii+1);
      wcs3D_[0] = wcs3D_[0] || (astGetI(ff, "Naxes")>2) ? 1 : 0;
    }
  }
  else {
    for (int ii=0; ii<nn; ii++) {
      AstFrame* ff = (AstFrame*)astGetFrame(ast_,ii+1);
      const char* id = astGetC(ff,"Ident");
      if (id && *id) {
	int jj = (*id == ' ') ? 0 : *id-'@';
	wcs3D_[jj] = (astGetI(ff, "Naxes")>2) ? 1 : 0;
      }
    }
  }

  astEnd;
}

void FitsImage::wcsHPXInit()
{
  wcsHPX_ =0;

  if (!ast_)
    return;

  char key[] = "CTYPE1 ";
  if (image_) {
    const char* str = image_->getKeyword(key);
    if (str)
      if (!strncmp(str+5,"HPX",3))
	wcsHPX_ =1;
  }
}

void FitsImage::wcsSizeInit()
{
  // init wcsSize_ array
  if (wcsSize_)
    delete [] wcsSize_;
  wcsSize_ =NULL;

  wcsSize_ = new double[MULTWCS];
  for (int ii=0; ii<MULTWCS; ii++)
    wcsSize_[ii] = calcWCSSize((Coord::CoordSystem)(ii+Coord::WCS));
}

int FitsImage::checkWCS(Vector& vv)
{
  // check for reasonable values
  return (fabs(vv[0]) < FLT_MAX &&
	  fabs(vv[1]) < FLT_MAX) ? 1 : 0;
}

int FitsImage::checkWCS(Vector3d& vv)
{
  // check for reasonable values
  return (fabs(vv[0]) < FLT_MAX &&
	  fabs(vv[1]) < FLT_MAX &&
	  fabs(vv[2]) < FLT_MAX ) ? 1 : 0;
}

void FitsImage::setWCSFormat(int id, const char* format)
{
  // is it already set?
  // ast is very slow when changing params
  {
    ostringstream str;
    str << "Format(" << id << ")" << ends;
    const char* out = astGetC(ast_, str.str().c_str());
    if (!strcmp(out,format))
      return;
  }

  ostringstream str;
  str << "Format(" << id << ")=" << format << ends;
  astSet(ast_, str.str().c_str());
}

void FitsImage::setWCSSkyFrame(Coord::CoordSystem sys, Coord::SkyFrame sky)
{
  int nn = astGetI(ast_,"nframe");
  char cc = ' ';
  int ww = sys-Coord::WCS;
  switch (sys) {
  case Coord::DATA:
  case Coord::IMAGE:
  case Coord::PHYSICAL:
  case Coord::AMPLIFIER:
  case Coord::DETECTOR:
    // this should not happen
    return;
  default:
    if (ww)
      cc = ww+'@';
  }

  int found=0;
  for (int ss=0; ss<nn; ss++) {
    const char* id = astGetC(astGetFrame(ast_,ss+1),"Ident");
    if (cc == id[0]) {
      astSetI(ast_,"Current",ss+1);
      found =1;
      break;
    }
  }

  // did we find it?
  if (!found)
    return; 
  
  // is celestial system?
  if (!hasWCSCel(sys))
    return;
  
  const char* str = astGetC(ast_, "System");
  switch (sky) {
  case Coord::FK4_NO_E:
    if (!strncmp(str,"FK4-NO-E",8))
      return;
    astSet(ast_, "System=FK4-NO-E, Equinox=B1950");
    return;
  case Coord::FK4:
    if (!strncmp(str,"FK4",3))
      return;
    astSet(ast_, "System=FK4, Equinox=B1950");
    return;
  case Coord::FK5:
    if (!strncmp(str,"FK5",3))
      return;
    astSet(ast_, "System=FK5, Equinox=J2000");
    return;
  case Coord::ICRS:
    if (!strncmp(str,"ICRS",4))
      return;
    astSet(ast_, "System=ICRS");
    return;
  case Coord::GALACTIC:
    if (!strncmp(str,"GALACTIC",8))
      return;
    astSet(ast_, "System=GALACTIC");
    return;
  case Coord::SUPERGALACTIC:
    if (!strncmp(str,"SUPERGALACTIC",13))
      return;
    astSet(ast_, "System=SUPERGALACTIC");
    return;
  case Coord::ECLIPTIC:
    if (!strncmp(str,"ECLIPTIC",8))
      return;
    astSet(ast_, "System=ECLIPTIC");
    // get AST to agree with WCSSUBS
    astSetD(ast_, "EQUINOX", astGetD(ast_, "EPOCH"));
    return;
  case Coord::HELIOECLIPTIC:
    if (!strncmp(str,"HELIOECLIPTIC",13))
      return;
    astSet(ast_, "System=HELIOECLIPTIC");
    return;
  }
}

Vector FitsImage::wcsTran(const Vector& in, int forward)
{
  int naxes = astGetI(ast_,"Naxes");
  switch (naxes) {
  case 1:
    // error
    break;
  case 2:
    double xout, yout;
    astTran2(ast_, 1, in.v, in.v+1, forward, &xout, &yout);
    return Vector(xout, yout);
  case 3:
    {
      double pin[3];
      double pout[3];
      pin[0] = in[0];
      pin[1] = in[1];
      pin[2] = forward ? context_->slice(2) : 0;
      astTranN(ast_, 1, 3, 1, pin, forward, 3, 1, pout);
      return Vector(pout[0],pout[1]);
    }
    break;
  case 4:
    {
      double pin[4];
      double pout[4];
      pin[0] = in[0];
      pin[1] = in[1];
      pin[2] = forward ? context_->slice(2) : 0;
      pin[3] = forward ? context_->slice(3) : 0;
      astTranN(ast_, 1, 4, 1, pin, forward, 4, 1, pout);
      return Vector(pout[0],pout[1]);
    }
    break;
  }
  return Vector();
}

void FitsImage::wcsTran(AstFrameSet* ast, int npoint, 
			Vector* in, int forward, Vector* out)
{
  int naxes = astGetI(ast,"Naxes");
  switch (naxes) {
  case 1:
    // error
    break;
  case 2:
    {
      double* xin = new double[npoint];
      double* yin = new double[npoint];
      double* xout = new double[npoint];
      double* yout = new double[npoint];

      for (int ii=0; ii<npoint; ii++) {
	xin[ii] = in[ii][0];
	yin[ii] = in[ii][1];
      }
      astTran2(ast, npoint, xin, yin, forward, xout, yout);
      for (int ii=0; ii<npoint; ii++)
	out[ii] = Vector(xout[ii],yout[ii]);

      if (xin)
	delete [] xin;
      if (yin)
	delete [] yin;
      if (xout)
	delete [] xout;
      if (yout)
	delete [] yout;
    }
    break;
  case 3:
    {
      double* ptr_in[3];
      ptr_in[0] = new double[npoint];
      ptr_in[1] = new double[npoint];
      ptr_in[2] = new double[npoint];
      double* ptr_out[3];
      ptr_out[0] = new double[npoint];
      ptr_out[1] = new double[npoint];
      ptr_out[2] = new double[npoint];

      for (int kk=0; kk<npoint; kk++) {
	ptr_in[0][kk] = in[kk][0];
	ptr_in[1][kk] = in[kk][1];
	ptr_in[2][kk] = forward ? context_->slice(2) : 0;
      }      
      astTranP(ast, npoint, 3, (const double**)ptr_in, forward, 3, ptr_out);
      for (int kk=0; kk<npoint; kk++)
	out[kk] = Vector(ptr_out[0][kk], ptr_out[1][kk]);

      if (ptr_in[0])
	delete [] ptr_in[0];
      if (ptr_in[1])
	delete [] ptr_in[1];
      if (ptr_in[2])
	delete [] ptr_in[2];

      if (ptr_out[0])
	delete [] ptr_out[0];
      if (ptr_out[1])
	delete [] ptr_out[1];
      if (ptr_out[2])
	delete [] ptr_out[2];
    }
    break;
  case 4:
    {
      double* ptr_in[4];
      ptr_in[0] = new double[npoint];
      ptr_in[1] = new double[npoint];
      ptr_in[2] = new double[npoint];
      ptr_in[3] = new double[npoint];
      double* ptr_out[4];
      ptr_out[0] = new double[npoint];
      ptr_out[1] = new double[npoint];
      ptr_out[2] = new double[npoint];
      ptr_out[3] = new double[npoint];

      for (int kk=0; kk<npoint; kk++) {
	ptr_in[0][kk] = in[kk][0];
	ptr_in[1][kk] = in[kk][1];
	ptr_in[2][kk] = forward ? context_->slice(2) : 0;
	ptr_in[3][kk] = forward ? context_->slice(3) : 0;
      }      
      astTranP(ast, npoint, 4, (const double**)ptr_in, forward, 4, ptr_out);
      for (int kk=0; kk<npoint; kk++)
	out[kk] = Vector(ptr_out[0][kk], ptr_out[1][kk]);

      if (ptr_in[0])
	delete [] ptr_in[0];
      if (ptr_in[1])
	delete [] ptr_in[1];
      if (ptr_in[2])
	delete [] ptr_in[2];
      if (ptr_in[3])
	delete [] ptr_in[3];

      if (ptr_out[0])
	delete [] ptr_out[0];
      if (ptr_out[1])
	delete [] ptr_out[1];
      if (ptr_out[2])
	delete [] ptr_out[2];
      if (ptr_out[3])
	delete [] ptr_out[3];
    }
    break;
  }
}

Vector3d FitsImage::wcsTran(const Vector3d& in, int forward)
{
  int naxes = astGetI(ast_,"Naxes");
  switch (naxes) {
  case 1:
  case 2:
      double pin[2];
      double pout[2];
      pin[0] = in[0];
      pin[1] = in[1];
      astTranN(ast_, 1, 2, 1, pin, forward, 2, 1, pout);
      return Vector3d(pout[0],pout[1],forward ? 1 : 0);
    break;
  case 3:
    {
      double pin[3];
      double pout[3];
      pin[0] = in[0];
      pin[1] = in[1];
      pin[2] = in[2];
      astTranN(ast_, 1, 3, 1, pin, forward, 3, 1, pout);
      return Vector3d(pout[0],pout[1],pout[2]);
    }
    break;
  case 4:
    {
      double pin[4];
      double pout[4];
      pin[0] = in[0];
      pin[1] = in[1];
      pin[2] = in[2];
      pin[3] = forward ? context_->slice(3) : 0;
      astTranN(ast_, 1, 4, 1, pin, forward, 4, 1, pout);
      return Vector3d(pout[0],pout[1],pout[2]);
    }
    break;
  }
  return Vector3d();
}

double FitsImage::wcsDistance(const Vector& vv1, const Vector& vv2)
{
  int naxes = astGetI(ast_,"Naxes");
  switch (naxes) {
  case 1:
    // error
    break;
  case 2:
    return astDistance(ast_, vv1.v, vv2.v);
  case 3:
    {
      double ptr1[3];
      ptr1[0] = vv1[0];
      ptr1[1] = vv1[1];
      ptr1[2] = 0;
      double ptr2[3];
      ptr2[0] = vv2[0];
      ptr2[1] = vv2[1];
      ptr2[2] = 0;

      return astDistance(ast_, ptr1, ptr2);
    }
  case 4:
    {
      double ptr1[4];
      ptr1[0] = vv1[0];
      ptr1[1] = vv1[1];
      ptr1[2] = 0;
      ptr1[3] = 0;
      double ptr2[4];
      ptr2[0] = vv2[0];
      ptr2[1] = vv2[1];
      ptr2[2] = 0;
      ptr2[3] = 0;

      return astDistance(ast_, ptr1, ptr2);
    }
  }

  return 0;
}

double FitsImage::wcsAngle(const Vector& vv1, const Vector& vv2,
			   const Vector& vv3)
{
  int naxes = astGetI(ast_,"Naxes");
  switch (naxes) {
  case 1:
    // error
    break;
  case 2:
    return astAngle(ast_,vv1.v,vv2.v,vv3.v);
  case 3:
    {
      double ptr1[3];
      ptr1[0] = vv1[0];
      ptr1[1] = vv1[1];
      ptr1[2] = 0;
      double ptr2[3];
      ptr2[0] = vv2[0];
      ptr2[1] = vv2[1];
      ptr2[2] = 0;
      double ptr3[3];
      ptr3[0] = vv3[0];
      ptr3[1] = vv3[1];
      ptr3[2] = 0;

      return astAngle(ast_, ptr1, ptr2, ptr3);
    }
  case 4:
    {
      double ptr1[4];
      ptr1[0] = vv1[0];
      ptr1[1] = vv1[1];
      ptr1[2] = 0;
      ptr1[3] = 0;
      double ptr2[4];
      ptr2[0] = vv2[0];
      ptr2[1] = vv2[1];
      ptr2[2] = 0;
      ptr2[3] = 0;
      double ptr3[4];
      ptr3[0] = vv3[0];
      ptr3[1] = vv3[1];
      ptr3[2] = 0;
      ptr3[3] = 0;

      return astAngle(ast_, ptr1, ptr2, ptr3);
    }
  }

  return 0;
}

double FitsImage::wcsAxAngle(const Vector& vv1, const Vector& vv2)
{
  int naxes = astGetI(ast_,"Naxes");
  switch (naxes) {
  case 1:
    // error
    break;
  case 2:
    return astAxAngle(ast_, vv1.v, vv2.v, 2);
  case 3:
    {
      double ptr1[3];
      ptr1[0] = vv1[0];
      ptr1[1] = vv1[1];
      ptr1[2] = 0;
      double ptr2[3];
      ptr2[0] = vv2[0];
      ptr2[1] = vv2[1];
      ptr2[2] = 0;

      return astAxAngle(ast_, ptr1, ptr2, 2);
    }
  case 4:
    {
      double ptr1[4];
      ptr1[0] = vv1[0];
      ptr1[1] = vv1[1];
      ptr1[2] = 0;
      ptr1[3] = 0;
      double ptr2[4];
      ptr2[0] = vv2[0];
      ptr2[1] = vv2[1];
      ptr2[2] = 0;
      ptr2[3] = 0;

      return astAxAngle(ast_, ptr1, ptr2, 2);
    }
  }

  return 0;
}

static void fits2TAB(AstFitsChan* chan, const char* extname,
		     int extver, int extlevel, int* status)
{
  FitsFile* ptr = ((FitsImage*)astChannelData)->fitsFile();
  FitsFile* ext =NULL;
  bool first=true;
  bool found=false;
  while (!found) {
    ext = new FitsMosaicNextMMapIncr(ptr);
    if (!first)
      delete ptr;
    first =false;
    
    // EOF?
    if (!ext || !ext->isValid()) {
      if (ext)
	delete ext;
      
      *status = 0;
      return;
    }

    if (!ext->isBinTable())
      break;

    const char* name = ext->extname();
    int ver = ext->extver();
    int level = ext->extlevel();
    
    if (name) {
      if (!strncmp(extname,name,7) && extver==ver && extlevel==level) {
	found =true;
	break;
      }
    }

    ptr = ext;
  }  

  // ok, found it
  astClearStatus; // just to make sure
  astBegin; // start memory management

  FitsHead* hd = ext->head();
  FitsBinTableHDU* hdu = (FitsBinTableHDU*)hd->hdu();
  int cols = hdu->cols();
  int rows = hdu->rows();
  int rowlen = hdu->width();

  // create fitstable
  AstFitsChan* header = astFitsChan(NULL, NULL, "");
  char* cards = hd->cards();
  int ncards = hd->ncard();

  for (int ii=0; ii<ncards; ii++) {
    char buf[81];
    strncpy(buf,cards+(ii*80),80);
    buf[80] = '\0';

    astPutFits(header, buf, 0);
  }
  AstFitsTable* table = (AstFitsTable*)astFitsTable(header,"");

  for (int ii=0; ii<cols; ii++) {
    FitsBinColumn* col = (FitsBinColumn*)hdu->find(ii);
    int width = col->width();
    int repeat = col->repeat();

    char* ptr = (char*)ext->data();
    unsigned char* data = new unsigned char[width*rows];
    memset(data,0,width*rows);
    switch (col->type()) {
    case 'I':
      for (int ii=0; ii<rows; ii++, ptr+=rowlen)
	for (int jj=0; jj<repeat; jj++) {
	  short vv = col->value(ptr,jj);
	  memcpy(data+ii*width+jj*2,&vv,2);
	}
      break;
    case 'J':
      for (int ii=0; ii<rows; ii++, ptr+=rowlen)
	for (int jj=0; jj<repeat; jj++) {
	  int vv = col->value(ptr,jj);
	  memcpy(data+ii*width+jj*4,&vv,4);
	}
      break;
    case 'E':
      for (int ii=0; ii<rows; ii++, ptr+=rowlen)
	for (int jj=0; jj<repeat; jj++) {
	  float vv = col->value(ptr,jj);
	  memcpy(data+ii*width+jj*4,&vv,4);
	}
      break;
    case 'D':
      for (int ii=0; ii<rows; ii++, ptr+=rowlen)
	for (int jj=0; jj<repeat; jj++) {
	  double vv = col->value(ptr,jj);
	  memcpy(data+ii*width+jj*8,&vv,8);
	}
      break;
    }
    astPutColumnData(table, col->ttype(), 0, width*rows, data);
    if (data)
      delete [] data;
  }  

  astPutTable(chan, table, extname);

  astEnd; // now, clean up memory
  if (ext)
    delete ext;

  *status = 1;
}

AstFrameSet* FitsImage::fits2ast(FitsHead* hd) 
{
  // we may have an error, just reset
  astClearStatus;

  // new fitschan
  AstFitsChan* chan = astFitsChan(NULL, NULL, "");
  if (!astOK || chan == AST__NULL)
    return NULL;

  // enable -TAB
  astSetI(chan,"TabOK",1);
  astPutChannelData(chan, this);
  astTableSource(chan, fits2TAB);

  // no warning messages
  astClear(chan,"Warnings");

  // fill up chan
  char* cards =NULL;
  int ncards =0;

  if (hd) {
    cards = hd->cards();
    ncards = hd->ncard();
  }

  if (cards == NULL || ncards == 0)
    return NULL;

  for (int i=0; i<ncards; i++) {
    char buf[81];
    strncpy(buf,cards+(i*80),80);
    buf[80] = '\0';

    astPutFits(chan, buf, 0);
    // sometimes, we get a bad parse, just ignore
    if (!astOK)
      astClearStatus;
  }

  // we may have an error, just reset
  astClearStatus;
  astClear(chan, "Card");

  // parse header
  AstFrameSet* frameSet = (AstFrameSet*)astRead(chan);

  // do we have anything?
  if (!astOK || frameSet == AST__NULL || 
      strncmp(astGetC(frameSet,"Class"), "FrameSet", 8))
    return NULL;

  // warn if no inverse
  astInv_ = astGetI(frameSet, "TranInverse");
  if (!astInv_)
    internalError("Warning: the WCS has no defined inverse. Some functionality may not be available.");
  
  return frameSet;
}