M3_format_LevelS.c

/*******************************************************************************
*   M3:  M3_format_LevelS.c                                                    *
*                                                                              *
*   Version 2.0 April 2007                                                     *
*                                                                              *
*   Copyright (C) 2004  C.M. Cantalupo                                         *
*                                                                              *
*   This program is free software; you can redistribute it and/or modify       *
*   it under the terms of the GNU General Public License as published by       *
*   the Free Software Foundation; either version 2 of the License, or          *
*   (at your option) any later version.                                        *
*                                                                              *
*   This program is distributed in the hope that it will be useful,            *
*   but WITHOUT ANY WARRANTY; without even the implied warranty of             *
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
*   GNU General Public License for more details.                               *
*                                                                              *
*   You should have received a copy of the GNU General Public License          *
*   along with this program; if not, write to the Free Software                *
*   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA  *
*                                                                              *
*******************************************************************************/

#include "M3.h"
#include "M3_internal.h"
#include "M3_gcp.h"
#include "M3_system.h"
#include "M3_format_LevelS.h"
#include "fitsio.h"
#include "fftw3.h"
#include <ctype.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#include "chealpix.h"
#include "rotateCoords.h"

#define M3_CONTENT_LENGTH 256
#define M3_NOISE_GEN_CORRFACTOR 8
#define M3_NOISE_GEN_MINIMUM 32768

void M3_File_parseFormatString_LevelS( M3_File *theFile, M3_FormatStruct_LevelS *out )
{
  /* Converts the format string from theFile into the format structure including setting default values. */

  char *s1, *s2;
  char errorString[M3_CONTENT_LENGTH];
  char errorString2[M3_CONTENT_LENGTH];
  int n;
  static char lastFormatString[M3_CONTENT_LENGTH] = {0};
  static M3_FormatStruct_LevelS lastFormatStruct;

  /* Build in a last call cache, as we will often be parsing the same format strings over and over again.  */
  if( strcmp( theFile->format, lastFormatString ) == 0 )
  {
    *out = lastFormatStruct;
    return;
  }

  /* Set the default values.  */
  out->template = 0;
  out->pixelType = M3_DEFAULT_PIXEL_TYPE;
  out->psiPol = 0;
  out->rotateDirection[0] = '\0';
  out->pixelScheme = 'r';
  strcpy( out->rotateEquinox, "J2000");
  out->simmission2006 = 0;
  out->header2007 = 0;
  out->numPeriod = 0;
  out->detectorID[0] = '\0';
  if( theFile->fileType == M3_TOD_FILE_TYPE )
    out->firstSample = INT64_MAX;
  else
    out->firstSample = 0;
  out->lastSample = INT64_MAX;
  out->corLength = 65536;
  out->nside = 0;
  out->spinAxis = '\0';
  out->LFI = 0;

  /* LevelS_pixelType=I_nside=256_pixelScheme=ring_psiPol=0_template_rotateDirection=radec2gal_rotateEquinox=J2000_simmission2006_numPeriod=61_header2007*/

  /* Set error string.  */
  sprintf( errorString, "Error interpreting format string: \"%s\"", theFile->format );

  /* Advance to the first parameter */
  s1 = strchr( theFile->format, '_');

  /* Loop over the paramters in the format string */
  while( s1 )
  {
    /* Skip '_' character */
    s1++;

    /* Switch over all of the parameters */
    if( strncmp( s1, "pixelType", 9 ) == 0 )
    {
      s1 += 9;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      if( strncmp( s1, "I", 1 ) == 0 )
        out->pixelType = M3_CMB_I_PIXEL_TYPE;
      else if( strncmp( s1, "Q", 1 ) == 0 )
        out->pixelType = M3_CMB_Q_PIXEL_TYPE;
      else if( strncmp( s1, "U", 1 ) == 0 )
        out->pixelType = M3_CMB_U_PIXEL_TYPE;
      else 
        out->pixelType = M3_DEFAULT_PIXEL_TYPE;      
    }
    else if( strncmp( s1, "nside", 5) == 0 )
    {
      s1 += 5;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      out->nside = atoi( s1 );
    }
    else if( strncmp( s1, "pixelScheme", 11 ) == 0 )
    {
      s1 += 11;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      if( strncmp( s1, "ring", 4) == 0 )
        out->pixelScheme = 'r';
      else if( strncmp( s1, "nest", 4) == 0 )
        out->pixelScheme = 'n';
      else
        M3_ErrorCheck( -1, errorString, 0, M3_EFLAG_DEFAULT );
    }
    else if( strncmp( s1, "psiPol", 6) == 0 )
    {
      s1 += 6;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      out->psiPol = atof(s1);
    }
    else if( strncmp( s1, "template", 8) == 0 )
    {
      out->template = 1;
    }
    else if( strncmp( s1, "rotateDirection", 15 ) == 0 )
    {
      s1 += 15;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      s2 = strchr( s1, '_');
      if( s2 == NULL )
        n = strlen(s1);
      else
        n = (size_t)(s2)-(size_t)(s1);

      sprintf( errorString2, "Error interpreting format string: \"%s\" rotateDirection too long", theFile->format ); 
      M3_ErrorCheck( -1, errorString2, n < 128, M3_EFLAG_DEFAULT );

      strncpy( out->rotateDirection, s1, n);
      out->rotateDirection[n] = '\0';
    }
    else if( strncmp( s1, "rotateEquinox", 13 ) == 0 )
    {
      s1 += 13;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      s2 = strchr( s1, '_');
      if( s2 == NULL )
        n = strlen(s1);
      else
        n = (size_t)(s2)-(size_t)(s1);
      
      sprintf( errorString2, "Error interpreting format string: \"%s\" rotateEquinox too long", theFile->format ); 
      M3_ErrorCheck( -1, errorString2, n < 128, M3_EFLAG_DEFAULT );

      strncpy( out->rotateEquinox, s1, n);
      out->rotateEquinox[n] = '\0';
    }
    else if( strncmp( s1, "simmission2006", 14 ) == 0 )
    {
      out->simmission2006 = 1;
    }
    else if( strncmp( s1, "numPeriod", 9 ) == 0 )
    {
      s1 += 9;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      out->numPeriod = atoi(s1);
    }
    else if( strncmp( s1, "header2007", 10 ) == 0 )
    {
      out->header2007 = 1;
    }
    else if( strncmp( s1, "firstSample", 11 ) == 0 )
    {
      s1 += 11;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->firstSample = strtoll(s1, NULL, 10 );
    }
    else if( strncmp( s1, "lastSample", 10 ) == 0 )
    {
      s1 += 10;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->lastSample = strtoll(s1, NULL, 10 );
    }
    else if( strncmp( s1, "corLength", 9 ) == 0 )
    {
      s1 += 9;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->corLength = strtoll(s1, NULL, 10 );
    }
    else if( strncmp( s1, "detectorID", 10 ) == 0 )
    {
      s1 += 10;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      s2 = strchr( s1, '_');
      if( s2 == NULL )
        n = strlen(s1);
      else
        n = (size_t)(s2)-(size_t)(s1);
      sprintf( errorString2, "Error interpreting format string: \"%s\" detectorID too long", theFile->format ); 
      M3_ErrorCheck( -1, errorString2, n < 128, M3_EFLAG_DEFAULT );
      strncpy( out->detectorID, s1, n);
      out->detectorID[n] = '\0';
    }
    else if( strncmp( s1, "detID", 5 ) == 0 )
    {
      s1 += 5;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      s2 = strchr( s1, '_');
      if( s2 == NULL )
        n = strlen(s1);
      else
        n = (size_t)(s2)-(size_t)(s1);
      sprintf( errorString2, "Error interpreting format string: \"%s\" detectorID too long", theFile->format ); 
      M3_ErrorCheck( -1, errorString2, n < 128, M3_EFLAG_DEFAULT );
      strncpy( out->detectorID, s1, n);
      out->detectorID[n] = '\0';
    }
    else if( strncmp( s1, "spinAxis", 8 ) == 0 )
    {
      s1 += 8;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      if( *s1 == 'x' || *s1 == 'X' )
        out->spinAxis = 'x';
      if( *s1 == 'z' || *s1 == 'Z' )
        out->spinAxis = 'z';
    }
    else if( strncmp( s1, "LFI", 3 ) == 0 )
    {
      out->LFI = 1;
    }

    s1 = strchr( s1, '_');
  }
  
  /* Set last call cache */
  strcpy( lastFormatString, theFile->format );
  lastFormatStruct = *out;
}





void M3_File_readHeader_LevelS( M3_File *theFile, void *data )
{
  /* Simple switch over sub-functions for each file type */
  switch( theFile->fileType )
  {
    case M3_TOD_FILE_TYPE:
      M3_File_readHeader_LevelS_tod( theFile, (M3_TODheader *)data );
      break;
    case M3_POINTING_FILE_TYPE:
      M3_File_readHeader_LevelS_pointing( theFile, (M3_PointingHeader *)data );
      break;
    case M3_GCPOINTING_FILE_TYPE:
      M3_File_readHeader_LevelS_gcpointing( theFile, (M3_GCPointingHeader *)data );
      break;
    case M3_AUX_FILE_TYPE:
      M3_File_readHeader_LevelS_aux( theFile, (M3_AuxHeader *)data );
      break;
    case M3_NOISE_FILE_TYPE:
      M3_File_readHeader_LevelS_noise( theFile, (M3_NoiseHeader *)data );
      break;
    default:
      M3_ErrorCheck(-1, "M3 can only read LevelS files for tod, pointing, gcpointing, and aux file types", 0, M3_EFLAG_DEFAULT );
      break;
  }
}


void M3_File_readData_LevelS( M3_File *theFile, void *readArg, void *data )
{
  /* Simple switch over sub-functions for each file type */
  switch( theFile->fileType )
  {
    case M3_TOD_FILE_TYPE:
      M3_File_readData_LevelS_tod( theFile, (M3_Interval*)readArg, (double *)data );
      break;
    case M3_POINTING_FILE_TYPE:
      M3_File_readData_LevelS_pointing( theFile, (M3_Interval*)readArg, (M3_PointingEl *)data );
      break;
    case M3_GCPOINTING_FILE_TYPE:
      M3_File_readData_LevelS_gcpointing( theFile, (M3_Interval*)readArg, (double *)data );
      break;
    case M3_AUX_FILE_TYPE:
      M3_File_readData_LevelS_aux( theFile, (M3_planckAux *)data );
      break;
    case M3_NOISE_FILE_TYPE:
      M3_File_readData_LevelS_noise( theFile, (double *)data );
      break;
    default:
      M3_ErrorCheck(-1, "M3 can only read LevelS files for tod, pointing, gcpointing, and aux file types", 0, M3_EFLAG_DEFAULT );
      break;
  }
}



void M3_File_readHeader_LevelS_tod( M3_File *theFile, M3_TODheader *data )
{
  int64_t firstPeriod, lastPeriod, samplesPerPeriod, numSample, naxis1;
  fitsfile *thisFitsFile;
  int status = 0;
  M3_FormatStruct_LevelS formatStruct;
  char *s1;
  char strBuffer[M3_CONTENT_LENGTH];
  char errorString[M3_CONTENT_LENGTH];

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status);
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  /* First look to see if "firstSample" was specified in format string, this over rides other derivations */  
  fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
  sprintf( errorString, "Error going to second fits header in file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );
  
  fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
  sprintf( errorString, "Could not read NAXIS2 flag from second header unit in 2007 style header in file %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  numSample = strtoll(strBuffer, NULL, 10);

  if( formatStruct.firstSample != INT64_MAX )
  {
    data->firstSample = formatStruct.firstSample;
    data->lastSample = formatStruct.firstSample + numSample - 1;

    fits_close_file( thisFitsFile, &status);
    sprintf( errorString, "Error closing fits file named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

    return;
  }

  /* Next check to see if we can find "first_sample" in the second header data unit.  */
  fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
  sprintf( errorString, "Error going to second fits header in file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );
  fits_read_key_str(thisFitsFile, "first_sample", strBuffer, NULL, &status);
  if( status == 0 )
  {
    data->firstSample = strtoll(strBuffer, NULL, 10);
    
    fits_read_key_str(thisFitsFile, "last_sample", strBuffer, NULL, &status);
    if( status == 0 )
      data->lastSample = strtoll(strBuffer, NULL, 10);
    else
    {
      data->lastSample = data->firstSample + numSample - 1;
      status = 0;
    }
    fits_close_file( thisFitsFile, &status);
    sprintf( errorString, "Error closing fits file named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

    return;
  }
  status = 0;

  /* If we can't find out about samples, figure out pointing periods.  */
  if( formatStruct.header2007 )
  {
    fits_movabs_hdu( thisFitsFile, 1, NULL, &status );
    sprintf( errorString, "Could not go to first fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "first_pointing", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read first_pointing flag from 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    firstPeriod = atoi( strBuffer );

    fits_read_key_str( thisFitsFile, "last_pointing", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read last_pointing flag from 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    lastPeriod = atoi( strBuffer );

    fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
    sprintf( errorString, "Error going to second fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read NAXIS2 flag from second header unit in 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    numSample = strtoll(strBuffer, NULL, 10);
  }
  else
  {
    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Error in skipping first fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    numSample = strtoll(strBuffer, NULL, 10);

    /* Fix for Torsti's sorption cooler files. */
    fits_read_key_str(thisFitsFile, "NAXIS1", strBuffer, NULL, &status);
    sprintf( errorString, "Error reading NAXIS1 key from fits header of file %s", theFile->name);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    naxis1 = strtoll(strBuffer, NULL, 10 );

    if( naxis1 == 48000 )
      numSample *= naxis1/4;
    
    status = 0;
    fits_read_key_str(thisFitsFile, "FIRSTPNT", strBuffer, NULL, &status);
    
    if( status != 0 )
    {
      sprintf( errorString, "FIRSTPNT key not found, and format string does not specify numPeriod in file %s", theFile->name );
      M3_ErrorCheck( -1, errorString, formatStruct.numPeriod, M3_EFLAG_DEFAULT );
      s1 = strstr( theFile->name, "_tod_" );
      if( s1 && isdigit( s1[5] ) )
      {
        firstPeriod = atoi(s1+5) + 1;
        lastPeriod = firstPeriod + formatStruct.numPeriod-1;
      }
      else 
      {
        s1 = strstr( theFile->name, ".fits");
        if( s1 )
        {
          s1--;
          while( s1 >= theFile->name && isdigit(*s1) )
            s1--;
          s1++;
          if( isdigit(*s1))
            firstPeriod = atoi(s1) + 1;
        }
        else
        {
          sprintf( errorString, "FIRSTPNT key not found and could not find pointing period in file name: %s", theFile->name );
          M3_ErrorCheck( -1, errorString, 0, M3_EFLAG_DEFAULT);
        }
      }
      status = 0;
    }
    else
    {
      firstPeriod = strtoll(strBuffer, NULL, 10);
    }
    
    fits_read_key_str(thisFitsFile, "LASTPNT", strBuffer, NULL, &status);
    if( status != 0 )
    {
      sprintf( errorString, "LASTPNT key not found, and format string does not specify numPeriod in file %s", theFile->name );
      M3_ErrorCheck( -1, errorString, formatStruct.numPeriod, M3_EFLAG_DEFAULT );
      status = 0;
      lastPeriod = firstPeriod + formatStruct.numPeriod - 1;
    }
    else
    {
      lastPeriod = strtoll(strBuffer, NULL, 10);
    }    
  }

  fits_close_file( thisFitsFile, &status);
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );


  if( firstPeriod != 1 )
  {
    sprintf( errorString, "Number of samples is not divisable  by number of periods in file named %s", theFile->name );
    M3_ErrorCheck( -1, errorString, numSample%(lastPeriod - firstPeriod + 1) == 0, M3_EFLAG_DEFAULT );
    samplesPerPeriod = numSample/(lastPeriod - firstPeriod + 1);

    data->firstSample = samplesPerPeriod*(firstPeriod - 1);
    data->lastSample = samplesPerPeriod*lastPeriod - 1;
  }
  else
  {
    data->firstSample = 0;
    data->lastSample = numSample - 1;
  }
}  


void M3_File_readHeader_LevelS_pointing( M3_File *theFile, M3_PointingHeader *data )
{
  int64_t firstPeriod, lastPeriod, samplesPerPeriod, numSample, naxis1;
  fitsfile *thisFitsFile;
  int status = 0;
  M3_FormatStruct_LevelS formatStruct;
  char *s1;
  char strBuffer[M3_CONTENT_LENGTH];
  char errorString[M3_CONTENT_LENGTH];

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status);
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );
  
  if( formatStruct.firstSample != INT64_MAX )
  {
    fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
    sprintf( errorString, "Error going to second fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read NAXIS2 flag from second header unit in 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    numSample = strtoll(strBuffer, NULL, 10);
    data->firstSample = formatStruct.firstSample;
    data->lastSample = formatStruct.firstSample + numSample - 1;
    data->numPixelInClass = formatStruct.nside*formatStruct.nside*12;
    data->numNZ = 1;

    fits_close_file( thisFitsFile, &status);
    sprintf( errorString, "Error closing fits file named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

    return;
  }

  if( formatStruct.header2007 )
  {
    fits_movabs_hdu( thisFitsFile, 1, NULL, &status );
    sprintf( errorString, "Could not go to first fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "first_pointing", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read first_pointing flag from 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    firstPeriod = strtoll( strBuffer, NULL, 10 );

    fits_read_key_str( thisFitsFile, "last_pointing", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read last_pointing flag from 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    lastPeriod = strtoll( strBuffer, NULL, 10 );

    fits_movabs_hdu( thisFitsFile, 2, NULL, &status);
    sprintf( errorString, "Error going to second fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read NAXIS2 flag from second header unit in 2007 style header in file %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    numSample = strtoll(strBuffer, NULL, 10);
  }
  else
  {
    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Error in skipping first fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    numSample = strtoll(strBuffer, NULL, 10);

    fits_read_key_str(thisFitsFile, "FIRSTPNT", strBuffer, NULL, &status);
    if( status != 0 )
    {
      sprintf( errorString, "FIRSTPNT key not found, and format string does not specify numPeriod in file %s", theFile->name );
      M3_ErrorCheck( -1, errorString, formatStruct.numPeriod, M3_EFLAG_DEFAULT );

      s1 = strstr( theFile->name, "_pointing_" );
      if( s1 && isdigit( s1[5] ) )
      {
        firstPeriod = strtoll(s1+5, NULL, 10) + 1;
        lastPeriod = firstPeriod + formatStruct.numPeriod-1;
      }
      else 
      {
        s1 = strstr( theFile->name, ".fits");
        if( s1 )
        {
          s1--;
          while( s1 >= theFile->name && isdigit(*s1) )
            s1--;
          s1++;
          if( isdigit(*s1))
            firstPeriod = strtoll(s1, NULL, 10) + 1;
        }
        else
        {
          sprintf( errorString, "FIRSTPNT key not found and could not find pointing period in file name: %s", theFile->name );
          M3_ErrorCheck( -1, errorString, 0, M3_EFLAG_DEFAULT);
        }
      }
      status = 0;
    }
    else
    {
      firstPeriod = strtoll(strBuffer, NULL, 10);
    }

    fits_read_key_str(thisFitsFile, "LASTPNT", strBuffer, NULL, &status);
    if( status != 0 )
    {
      sprintf( errorString, "LASTPNT key not found, and format string does not specify numPeriod in file %s", theFile->name );
      M3_ErrorCheck( -1, errorString, formatStruct.numPeriod, M3_EFLAG_DEFAULT );
      status = 0;
      lastPeriod = firstPeriod + formatStruct.numPeriod - 1;
    }
    else
    {
      lastPeriod = strtoll(strBuffer, NULL, 10);
    }
  }

  fits_close_file( thisFitsFile, &status);
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

  if( firstPeriod != 1 )
  {
    sprintf( errorString, "Number of samples is not divisable  by number of periods in file named %s", theFile->name );
    M3_ErrorCheck( -1, errorString, numSample%(lastPeriod - firstPeriod + 1) == 0, M3_EFLAG_DEFAULT );
    samplesPerPeriod = numSample/(lastPeriod - firstPeriod + 1);

    data->firstSample = samplesPerPeriod*(firstPeriod - 1);
    data->lastSample = samplesPerPeriod*lastPeriod - 1;
  }
  else
  {
    data->firstSample = 0;
    data->lastSample = numSample - 1;
  }

  data->numPixelInClass = formatStruct.nside*formatStruct.nside*12;
  data->numNZ = 1;
}  

time_t my_timegm (struct tm *tm) 
{
  time_t ret;
  char *tz;
  
  tz = getenv("TZ");
  setenv("TZ", "", 1);
  tzset();
  ret = mktime(tm);
  if (tz)
    setenv("TZ", tz, 1);
  else
    unsetenv("TZ");
  tzset();
  return ret;
}

void M3_TimeEl_convertKeys( M3_TimeEl *outTime, char *dateKey, char *timeKey ) 
{
  char tempBuffer[M3_CONTENT_LENGTH];
  char *strPtr;
  struct tm timeStruct;

  strPtr = dateKey;

  memcpy( tempBuffer, strPtr, 4);
  tempBuffer[4] = 0;
  timeStruct.tm_year = atoi(tempBuffer) - 1900;
    
  memcpy( tempBuffer, strPtr+4, 2);
  tempBuffer[2] = 0;
  timeStruct.tm_mon = atoi(tempBuffer) - 1;
    
  memcpy( tempBuffer, strPtr+6, 2);
  tempBuffer[2] = 0;
  timeStruct.tm_mday = atoi(tempBuffer);
    
  strPtr = timeKey;

  memcpy( tempBuffer, strPtr, 2);
  tempBuffer[2] = 0;
  timeStruct.tm_hour = atoi( tempBuffer );
  
  memcpy( tempBuffer, strPtr+2, 2);
  tempBuffer[2] = 0;
  timeStruct.tm_min = atoi( tempBuffer );
  
  strcpy( tempBuffer, strPtr+4);
  timeStruct.tm_sec = atoi( tempBuffer );

  outTime->sec = my_timegm(&timeStruct);        
  outTime->nsec = 1e9*(atof(tempBuffer) - timeStruct.tm_sec);
}


void M3_File_readHeader_LevelS_gcpointing( M3_File *theFile, M3_GCPointingHeader *data )
{
  M3_FormatStruct_LevelS formatStruct;
  fitsfile *thisFitsFile;
  int status = 0;
  char errorString[M3_CONTENT_LENGTH];
  char temp[M3_CONTENT_LENGTH];
  char dateKey[M3_CONTENT_LENGTH];
  char timeKey[M3_CONTENT_LENGTH];
  char strBuffer[M3_CONTENT_LENGTH];
  char *strPtr;
  int64_t firstPeriod, lastPeriod;
  struct tm timeStruct;
  M3_TimeEl firstTime, secondTime, lastTime;
  int64_t numSample, numRow, n12;
  double tempd, sampleRate;

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status);
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  if( formatStruct.simmission2006 ) 
  {
    fits_movabs_hdu( thisFitsFile, 3, NULL, &status );
    sprintf( errorString, "Could not go to third fits HDU in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );
  } 
  else 
  {
    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Could not go to second fits HDU in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT ); 
  }
    
  fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
  sprintf( errorString, "Could not read NAXIS2 flag from third header unit in file:  %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  numSample = strtoll(strBuffer, NULL, 10);
  
  if( formatStruct.simmission2006 )
  {
    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Could not go to second fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str( thisFitsFile, "NAXIS2", strBuffer, NULL, &status );
    sprintf( errorString, "Could not read NAXIS2 flag from second header unit of file: %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    numRow = strtoll(strBuffer, NULL, 10);

    fits_read_col( thisFitsFile, TDOUBLE, 1, 1, 1, 1, NULL, &tempd, NULL, &status );
    firstTime.sec = tempd;
    fits_read_col( thisFitsFile, TDOUBLE, 2, 1, 1, 1, NULL, &tempd, NULL, &status );
    firstTime.nsec = tempd*1e9;

    fits_read_col( thisFitsFile, TDOUBLE, 13, 1, 1, 1, NULL, &tempd, NULL, &status );
    secondTime.sec = tempd;
    fits_read_col( thisFitsFile, TDOUBLE, 14, 1, 1, 1, NULL, &tempd, NULL, &status );
    secondTime.nsec = tempd*1e9;

    fits_read_col( thisFitsFile, TDOUBLE, 13, numRow, 1, 1, NULL, &tempd, NULL, &status );
    lastTime.sec = tempd;
    fits_read_col( thisFitsFile, TDOUBLE, 14, numRow, 1, 1, NULL, &tempd, NULL, &status );
    lastTime.nsec = tempd*1e9;

    /* If the file has variable pointing periods, then use the first period 
       to determine the rate to avoid using times with round off errors.  */
    fits_read_col( thisFitsFile, TLONGLONG, 18, 1, 1, 1, NULL, &n12, NULL, &status );
    if( status == 0 )
      sampleRate = n12/M3_TimeEl_Difference( secondTime, firstTime );
    else
    {
      status = 0;
      sampleRate = numSample/M3_TimeEl_Difference( lastTime, firstTime );
    }
  }
  else
  {
    fits_read_key_str(thisFitsFile, "FIRSTPT", strBuffer, NULL, &status);
    if( status != 0 )
    {
      firstPeriod = 0;
      status = 0;
    }
    fits_read_key_str(thisFitsFile, "LASTPT", strBuffer, NULL, &status);
    if( status != 0 )
    {
      lastPeriod = 8783;
      status = 0;
    }

    sprintf(temp, "SI%6lld", firstPeriod);
    fits_read_key_str(thisFitsFile, temp, strBuffer, NULL, &status );
    sprintf(errorString, "Error reading %s key from fits header of file %s", temp, theFile->name);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
    firstTime.sec = (int64_t)atof(strBuffer);

    sprintf(temp, "SF%6lld", firstPeriod);
    fits_read_key_str(thisFitsFile, temp, strBuffer, NULL, &status );
    sprintf(errorString, "Error reading %s key from fits header of file %s", theFile->name, temp);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);      
    firstTime.nsec = 1e9*atof(strBuffer);

    sprintf(temp, "EI%6lld", lastPeriod);
    fits_read_key_str(thisFitsFile, temp, strBuffer, NULL, &status );
    sprintf(errorString, "Error reading %s key from fits header of file %s", temp, theFile->name);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
    lastTime.sec = (int64_t)atof(strBuffer);
    
    sprintf(temp, "EF%6lld", lastPeriod);
    fits_read_key_str(thisFitsFile, temp, strBuffer, NULL, &status );
    sprintf(errorString, "Error reading %s key from fits header of file %s", temp, theFile->name);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT);
    lastTime.nsec = 1e9*atof(strBuffer);

    sampleRate = numSample/M3_TimeEl_Difference( lastTime, firstTime );

  }

  fits_close_file( thisFitsFile, &status);
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

  data->firstTime_sec = firstTime.sec;
  data->firstTime_nsec = firstTime.nsec;
  data->numSample = numSample;
  data->sampleRate = sampleRate;
  data->numDataPerSample = 6;

}

void M3_File_readHeader_LevelS_aux( M3_File *theFile, M3_AuxHeader *data )
{
  int64_t tableLength, tableWidth;
  M3_FormatStruct_LevelS formatStruct;
  fitsfile *thisFitsFile;
  int status = 0;
  char errorString[M3_CONTENT_LENGTH];
  char strBuffer[M3_CONTENT_LENGTH];

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status);
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
  sprintf( errorString, "Could not go to second fits header in file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  fits_read_key_str(thisFitsFile, "NAXIS1", strBuffer, NULL, &status);
  sprintf( errorString, "Error reading NAXIS1 key from second HDU of fits header of file %s", theFile->name);  
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  tableWidth = strtoll(strBuffer, NULL, 10);

  fits_read_key_str(thisFitsFile, "NAXIS2", strBuffer, NULL, &status);
  sprintf( errorString, "Error reading NAXIS2 key from second HDU of fits header of file %s", theFile->name);
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  tableLength = strtoll(strBuffer, NULL, 10);

  data->size = sizeof(M3_planckAux) + tableLength*(sizeof(M3_planckDetector) + 8);

  fits_close_file( thisFitsFile, &status);
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
}


void M3_File_readHeader_LevelS_noise( M3_File *theFile, M3_NoiseHeader *header )
{
  M3_FormatStruct_LevelS formatStruct;

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );
  header->firstSample = formatStruct.firstSample;
  header->lastSample = formatStruct.lastSample;
  header->corLength = formatStruct.corLength;
}

void M3_File_readData_LevelS_tod( M3_File *theFile, M3_Interval *readInterval, double *data )
{
  M3_FormatStruct_LevelS formatStruct;
  float *fltPtr;
  int64_t readLength, offset, i;
  FILE *fid;
  int isBigEndian;
  int status = 0;
  fitsfile *thisFitsFile;
  int colnums;
  char errorString[M3_CONTENT_LENGTH];

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  if (formatStruct.LFI) {

    colnums = 3;

    fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );

    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Could not go to second fits HDU in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    status = read_fits_tableblock(thisFitsFile, readInterval->firstSample, readInterval->lastSample, NULL, 1, FITS_BLOCK_NUMS, NULL, &colnums, data );

    sprintf( errorString, "Could not read samples %lld to %lld of LevelS TOD %s", readInterval->firstSample, readInterval->lastSample, theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_close_file( thisFitsFile, &status );
    sprintf( errorString, "Error closing fits file named %s", theFile->name );
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

  } else {

    fid = fopen( theFile->name, "r");
    M3_ErrorCheck( -1, theFile->name, fid, M3_EFLAG_FOPEN_READ );
    skipFitsHeaders( fid, 2);
    isBigEndian = testEndian();

    readLength = readInterval->lastSample - readInterval->firstSample + 1;
    offset = sizeof(float)*(readInterval->firstSample - theFile->param.tod.interval.firstSample );
    fseeko64(fid, offset, SEEK_CUR);
    fltPtr = ((float*)data) + readLength;
    fread(fltPtr, sizeof(float), readLength, fid);
    if( !isBigEndian )
      flipByteOrder( fltPtr, readLength, sizeof(float));
    for( i = 0; i < readLength; i++ )
      data[i] = (double)(fltPtr[i]);
    fclose(fid);

  }

  return;
}


#ifdef UNDEFINED
void M3_File_readData_LevelS_tod( M3_File *theFile, M3_Interval *readInterval, double *data )
{
  float *fltPtr;
  int64_t readLength, offset, i;
  FILE *fid;
  int isBigEndian;
  int status = 0;
  fitsfile *thisFitsFile;
  int colnums = 0;
  char errorString[M3_CONTENT_LENGTH];
  M3_FormatStruct_LevelS formatStruct;

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );

  fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
  sprintf( errorString, "Could not go to second fits HDU in file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );
    
  if( formatStruct.LFI )
    colnums = 3;
  else
    colnums = 0;

  status = read_fits_tableblock(thisFitsFile, readInterval->firstSample, readInterval->lastSample, NULL, 1, FITS_BLOCK_NUMS, NULL, &colnums, data );

  sprintf( errorString, "Could not read samples %lld to %lld of LevelS TOD %s", readInterval->firstSample, readInterval->lastSample, theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  fits_close_file( thisFitsFile, &status );
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  return;
}
#endif

void M3_File_readData_LevelS_pointing( M3_File *theFile, M3_Interval *readInterval,  M3_PointingEl *data )
{
  M3_FormatStruct_LevelS formatStruct;
  int64_t readLength, offset, i;
  FILE *fid;
  fitsfile *thisFitsFile;
  int status = 0;
  char errorString[M3_CONTENT_LENGTH];
  char strBuffer[M3_CONTENT_LENGTH];
  int coordBytes, isBigEndian;
  int64_t bufferLength, numReads, j;
  double coordDblPtr[3];
  void *fileBuffer;
  float *fltPtr;
  double *dblPtr;
  M3_PointingEl *pointingPtr;
  long p, nside;

  static int64_t diskBufferSize = -1;
  char *tempPtr;

  if( diskBufferSize == -1 )
  {
    tempPtr = getenv("M3_DISK_BUFFER_SIZE");
    if( tempPtr )
      diskBufferSize = strtoll( tempPtr, NULL, 10);
    else
      diskBufferSize = M3_DISK_BUFFER_SIZE;
  }



  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  /* Use cfitsio to determine if the data is double or single precision */
  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
  sprintf( errorString, "Could not go to second fits header in file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  fits_read_key_str( thisFitsFile, "TFORM1", strBuffer, NULL, &status );
  sprintf( errorString, "Could not read TFORM1 flag from second HDU in file %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  if( strncmp( strBuffer, "1E", 2) == 0  || strncmp( strBuffer+1, "1E", 2) == 0 )
    coordBytes = 4;
    else if( strncmp( strBuffer, "1D", 2) == 0  || strncmp( strBuffer+1, "1D", 2) == 0 )
    coordBytes = 8;
  else
  {
    sprintf( errorString, "Could not interpret TFORM1 flag from second HDU in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, 0, M3_EFLAG_DEFAULT );
  }

  fits_close_file( thisFitsFile, &status);
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

  fid = fopen( theFile->name, "r");
  M3_ErrorCheck( -1, theFile->name, fid, M3_EFLAG_FOPEN_READ );
  skipFitsHeaders( fid, 2);
  isBigEndian = testEndian();

  if( formatStruct.template )
  {
    readLength = readInterval->lastSample - readInterval->firstSample + 1;
    offset = sizeof(float)*(readInterval->firstSample - theFile->param.pointing.interval.firstSample );
    fseeko64( fid, offset, SEEK_CUR);
        
    fltPtr = ((float*)data) + readLength;
    fread(fltPtr, sizeof(float), readLength, fid);
    if( !isBigEndian )
      flipByteOrder( fltPtr, readLength, sizeof(float));

    for( i = 0; i < readLength; i++ )
    {
      data[i].pixel = 0;
      data[i].weight = fltPtr[i];
    }
  }
  else
  {
    fileBuffer = malloc(diskBufferSize);
    M3_ErrorCheck( -1, "M3_File_readData_LevelS_pointing()", (fileBuffer ? 1 : 0), M3_EFLAG_MALLOC_FSTRING );

    readLength = readInterval->lastSample - readInterval->firstSample + 1;
    bufferLength = diskBufferSize/(coordBytes*3);
    numReads = readLength/bufferLength + ( readLength%bufferLength ? 1 : 0 );

    offset = 3*coordBytes*(readInterval->firstSample - theFile->param.pointing.interval.firstSample );
    fseeko64( fid, offset, SEEK_CUR);

    fltPtr = (float *)fileBuffer;
    dblPtr = (double *)fileBuffer;
    pointingPtr = data;

    for( i = 0; i < numReads; i++ )
    {
      if( i == numReads - 1 && readLength%bufferLength )
        bufferLength = readLength%bufferLength;
      
      fread(fileBuffer, coordBytes, 3*bufferLength, fid );
      if( !isBigEndian )
        flipByteOrder( fileBuffer, 3*bufferLength, coordBytes );
      if( formatStruct.rotateDirection[0] != '\0' )
        if( coordBytes == 4 )
          for( j = 0; j < bufferLength; j++ )
          {
            coordDblPtr[0] = fltPtr[3*j];
            coordDblPtr[1] = fltPtr[3*j+1];
            coordDblPtr[2] = fltPtr[3*j+2];
            rotateCoords( formatStruct.rotateDirection, formatStruct.rotateEquinox, 1, coordDblPtr );
            fltPtr[3*j] = coordDblPtr[0];
            fltPtr[3*j+1] = coordDblPtr[1];
            fltPtr[3*j+2] = coordDblPtr[2];
          }
        else
          rotateCoords( formatStruct.rotateDirection, formatStruct.rotateEquinox, bufferLength, dblPtr );

      if( coordBytes  == 4 )
      {
        if( formatStruct.pixelScheme == 'r' )
        {
          if( formatStruct.pixelType == M3_CMB_I_PIXEL_TYPE || formatStruct.pixelType == M3_DEFAULT_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_ring(formatStruct.nside, fltPtr[3*j], fltPtr[3*j+1], &p);
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = 1;
            }
          else if( formatStruct.pixelType == M3_CMB_Q_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_ring(formatStruct.nside, fltPtr[3*j], fltPtr[3*j+1], &p );
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = cos(2*(fltPtr[3*j+2] + formatStruct.psiPol));
            }
          else if( formatStruct.pixelType == M3_CMB_U_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_ring(formatStruct.nside, fltPtr[3*j], fltPtr[3*j+1], &p );
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = sin(2*(fltPtr[3*j+2] + formatStruct.psiPol));
            }
        }
        else
        {
          if( formatStruct.pixelType == M3_CMB_I_PIXEL_TYPE || formatStruct.pixelType == M3_DEFAULT_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_nest(formatStruct.nside, fltPtr[3*j], fltPtr[3*j+1], &p);
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = 1;
            }
          else if( formatStruct.pixelType == M3_CMB_Q_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_nest(formatStruct.nside, fltPtr[3*j], fltPtr[3*j+1], &p );
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = cos(2*(fltPtr[3*j+2] + formatStruct.psiPol));
            }
          else if( formatStruct.pixelType == M3_CMB_U_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_nest(formatStruct.nside, fltPtr[3*j], fltPtr[3*j+1], &p );
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = sin(2*(fltPtr[3*j+2] + formatStruct.psiPol));
            }
        }
      }
      else 
      {
        if( formatStruct.pixelScheme == 'r' )
        {
          if( formatStruct.pixelType == M3_CMB_I_PIXEL_TYPE || formatStruct.pixelType == M3_DEFAULT_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_ring(formatStruct.nside, dblPtr[3*j], dblPtr[3*j+1], &p);
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = 1;
            }
          else if( formatStruct.pixelType == M3_CMB_Q_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_ring(formatStruct.nside, dblPtr[3*j], dblPtr[3*j+1], &p );
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = cos(2*(dblPtr[3*j+2] + formatStruct.psiPol));
            }
          else if( formatStruct.pixelType == M3_CMB_U_PIXEL_TYPE )
            for( j = 0; j < bufferLength; j++ )
            {
              ang2pix_ring(formatStruct.nside, dblPtr[3*j], dblPtr[3*j+1], &p );
              pointingPtr[j].pixel = p;
              pointingPtr[j].weight = sin(2*(dblPtr[3*j+2] + formatStruct.psiPol));
            }
          }
          else
          {
            if( formatStruct.pixelType == M3_CMB_I_PIXEL_TYPE || formatStruct.pixelType == M3_DEFAULT_PIXEL_TYPE )
              for( j = 0; j < bufferLength; j++ )
              {
                ang2pix_nest(formatStruct.nside, dblPtr[3*j], dblPtr[3*j+1], &p);
                pointingPtr[j].pixel = p;
                pointingPtr[j].weight = 1;
              }
            else if( formatStruct.pixelType == M3_CMB_Q_PIXEL_TYPE )
              for( j = 0; j < bufferLength; j++ )
              {
                ang2pix_nest(formatStruct.nside, dblPtr[3*j], dblPtr[3*j+1], &p );
                pointingPtr[j].pixel = p;
                pointingPtr[j].weight = cos(2*(dblPtr[3*j+2] + formatStruct.psiPol));
              }
            else if( formatStruct.pixelType == M3_CMB_U_PIXEL_TYPE )
              for( j = 0; j < bufferLength; j++ )
              {
                ang2pix_nest(formatStruct.nside, dblPtr[3*j], dblPtr[3*j+1], &p );
                pointingPtr[j].pixel = p;
                pointingPtr[j].weight = sin(2*(dblPtr[3*j+2] + formatStruct.psiPol));
              }
          }
      }
        pointingPtr += bufferLength;
    }
    free( fileBuffer );
  }
  fclose(fid);
}


void M3_File_readData_LevelS_gcpointing( M3_File *theFile, M3_Interval *readInterval, double *data )
{
  int64_t readLength;
  M3_FormatStruct_LevelS formatStruct;
  int64_t offset, i, j, k, l;
  double tempd;
  fitsfile *theFitsFile;
  int status = 0;
  char errorString[M3_CONTENT_LENGTH];
  char *tempPtr;
  int colnums[6] = {0, 1, 2, 3, 4, 5};


  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  readLength = readInterval->lastSample - readInterval->firstSample + 1;

  fits_open_file( &theFitsFile, theFile->name, READONLY, &status );
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );
  
  if( formatStruct.simmission2006 ) {

    fits_movabs_hdu( theFitsFile, 3, NULL, &status );
    sprintf( errorString, "Could not go to third fits HDU in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  } else {
    
    fits_movabs_hdu( theFitsFile, 2, NULL, &status );
    sprintf( errorString, "Could not go to second fits HDU in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );
    
  }
    
  status = read_fits_tableblock(theFitsFile, readInterval->firstSample, readInterval->lastSample, NULL, 6, FITS_BLOCK_NUMS, NULL, colnums, data);
    
  sprintf( errorString, "Error with block-reading file %s", theFile->name );
  M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  if( formatStruct.spinAxis != 'z' ) {
    for( i = 0; i < readLength; i++ ) {
      offset = 6*i;
      tempd = data[offset];
      data[offset] = data[offset+4];
      data[offset+4] = tempd;
      tempd = data[offset+1];
      data[offset+1] = data[offset+5];
      data[offset+5] = tempd;
      data[offset+2] = M_PI - data[offset+2];
      data[offset+3] = data[offset+3] + M_PI;
    }
  }
    
  fits_close_file(theFitsFile, &status);
  
  return;
}


void M3_File_readData_LevelS_aux( M3_File *theFile, M3_planckAux *data )
{
  void *fpBuffer;
  M3_FormatStruct_LevelS formatStruct;
  fitsfile *thisFitsFile;
  int status = 0, isBigEndian;
  char errorString[M3_CONTENT_LENGTH];
  char strBuffer[M3_CONTENT_LENGTH];
  double boresightSpinAxisAngle;  
  int64_t tableWidth, tableLength, i, j;
  FILE *fid;
  double *dblPtr;


  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  if( formatStruct.header2007 )
  {
    fits_movabs_hdu( thisFitsFile, 1, NULL, &status );
    sprintf( errorString, "Could not go to first fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str(thisFitsFile, "theta_b", strBuffer, NULL, &status);
    sprintf( errorString, "Error reading theta_b key from fits header of file %s", theFile->name);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    boresightSpinAxisAngle = (M_PI/180.0)*atof(strBuffer);

    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Could not go to second fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

  }
  else
  {
    fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
    sprintf( errorString, "Could not go to second fits header in file %s", theFile->name );
    M3_ErrorCheck( -1, errorString, status == 0, M3_EFLAG_DEFAULT );

    fits_read_key_str(thisFitsFile, "THETA_B", strBuffer, NULL, &status);
    sprintf( errorString, "Error reading THETA_B key from fits header of file %s", theFile->name);
    M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
    boresightSpinAxisAngle = (M_PI/180.0)*atof(strBuffer);
  }
  
  fits_read_key_str(thisFitsFile, "NAXIS1", strBuffer, NULL, &status);
  sprintf( errorString, "Error reading NAXIS1 key from second HDU of fits header of file %s", theFile->name);  
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  tableWidth = strtoll(strBuffer, NULL, 10);

  fits_read_key_str(thisFitsFile, "NAXIS2", strBuffer, NULL, &status);
  sprintf( errorString, "Error reading NAXIS2 key from second HDU of fits header of file %s", theFile->name);
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );
  tableLength = strtoll(strBuffer, NULL, 10);

  fits_close_file( thisFitsFile, &status);
  sprintf( errorString, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorString, status == 0, M3_EFLAG_DEFAULT );

  data->boresightSpinAxisAngle = boresightSpinAxisAngle;
  data->numDetector = tableLength;
  data->detectorInfo = (M3_planckDetector*)(((char*)data) + sizeof(M3_planckAux));


  fid = fopen( theFile->name, "r");
  M3_ErrorCheck( -1, theFile->name, fid, M3_EFLAG_FOPEN_READ );
  skipFitsHeaders( fid, 2);
  isBigEndian = testEndian();

  fpBuffer = malloc( tableWidth*tableLength);
  M3_ErrorCheck(-1, "M3_File_readData_LevelS_aux", (fpBuffer ? 1 : 0), M3_EFLAG_MALLOC_FSTRING);
  fread(fpBuffer, 1, tableWidth*tableLength, fid );
      
  for( i = 0; i < tableLength; i++ )
  {
    ((M3_planckAux*)data)->detectorInfo[i].detectorID = ((char*)data) + sizeof(M3_planckAux) + 
                                                        tableLength*sizeof(M3_planckDetector) + 
                                                        8*i;
    memcpy(((M3_planckAux*)data)->detectorInfo[i].detectorID, ((char*)fpBuffer) + i*tableWidth, 8);
    for( j = 0; j < 8; j++ )
      if( ((M3_planckAux*)data)->detectorInfo[i].detectorID[j] == ' ')
        break;
        
    ((M3_planckAux*)data)->detectorInfo[i].detectorID[j] = '\0';
        
    dblPtr = (double*)(((char*)fpBuffer) + i*tableWidth + 8);
    if(!isBigEndian)
      flipByteOrder( dblPtr, 4, sizeof(double));
    data->detectorInfo[i].offset[0] = (M_PI/180.0)*dblPtr[1];
    data->detectorInfo[i].offset[1] = (M_PI/180.0)*dblPtr[0];
    data->detectorInfo[i].offset[2] = (M_PI/180.0)*dblPtr[2];
    data->detectorInfo[i].psiPol = (M_PI/180.0)*dblPtr[3];
  }
  free(fpBuffer);

  fclose(fid);
}



void M3_File_getFPDBparam( M3_File *theFile, double *fknee, double *alpha, double *fmin, double *fsample, double *NET )
{
  M3_FormatStruct_LevelS formatStruct;
  fitsfile *thisFitsFile;
  int status = 0;
  long nrows;
  int colnum;
  char **detids;
  int coltype;
  long colrepeat;
  long colwidth;
  long detrow;

  int i, j, k, foundDet = -1;

  char errorMessage[M3_CONTENT_LENGTH];

  M3_File_parseFormatString_LevelS( theFile, &formatStruct );

  fits_open_file( &thisFitsFile, theFile->name, READONLY, &status );
  M3_ErrorCheck(-1, theFile->name, status == 0, M3_EFLAG_FOPEN_READ );

  fits_movabs_hdu( thisFitsFile, 2, NULL, &status );
  sprintf( errorMessage, "Could not go to second HDU in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  /* get the number of rows in the table */

  fits_get_num_rows( thisFitsFile, &nrows, &status);
  sprintf( errorMessage, "Could not retrieve number of rows in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  /* determine detector column and type */

  fits_get_colnum( thisFitsFile, CASEINSEN, "DETECTOR", &colnum, &status);
  sprintf( errorMessage, "Could not find detector column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  fits_get_coltype( thisFitsFile, colnum, &coltype, &colrepeat, &colwidth, &status);
  sprintf( errorMessage, "Could not determine type of detector column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  /* allocate memory for detector column */

  detids = (char**)malloc(nrows * sizeof(char*));
  sprintf( errorMessage, "Could not allocate detector IDs for file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  for (i = 0; i < nrows; i++) {
    detids[i] = (char*)malloc((colwidth+2) * sizeof(char));
    sprintf( errorMessage, "Could not allocate detector ID for file %s, row %d", theFile->name, i );
    M3_ErrorCheck( -1, errorMessage, detids[i], M3_EFLAG_DEFAULT );
  }
  
  /* read detector IDs and search for the row we want */

  fits_read_col( thisFitsFile, TSTRING, colnum, 1, 1, nrows, NULL, detids, NULL, &status);
  sprintf( errorMessage, "Could not read detector IDs for file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  for (i = 0; i < nrows; i++) {
    if( strcmp(detids[i], formatStruct.detectorID) == 0 ) {
      foundDet = i;
      break;
    }
  }
  if (foundDet == -1) {
    sprintf( errorMessage, "Could not find detector ID %s in file %s", formatStruct.detectorID, theFile->name );
    M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );
  }

  /* free detector ID memory */

  for (i = 0; i < nrows; i++) {
    free(detids[i]);
  }
  free(detids);

  /* for each parameter we want, look up column number and read the parameter */

  fits_get_colnum( thisFitsFile, CASEINSEN, "F_KNEE", &colnum, &status);
  sprintf( errorMessage, "Could not find f_knee column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  fits_read_col( thisFitsFile, TDOUBLE, colnum, (foundDet+1), 1, 1, NULL, fknee, NULL, &status);
  sprintf( errorMessage, "Could not read f_knee for detector %s in file %s", formatStruct.detectorID, theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  fits_get_colnum( thisFitsFile, CASEINSEN, "ALPHA", &colnum, &status);
  sprintf( errorMessage, "Could not find alpha column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  fits_read_col( thisFitsFile, TDOUBLE, colnum, (foundDet+1), 1, 1, NULL, alpha, NULL, &status);
  sprintf( errorMessage, "Could not read alpha for detector %s in file %s", formatStruct.detectorID, theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  fits_get_colnum( thisFitsFile, CASEINSEN, "F_MIN", &colnum, &status);
  sprintf( errorMessage, "Could not find f_min column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  fits_read_col( thisFitsFile, TDOUBLE, colnum, (foundDet+1), 1, 1, NULL, fmin, NULL, &status);
  sprintf( errorMessage, "Could not read f_min for detector %s in file %s", formatStruct.detectorID, theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  fits_get_colnum( thisFitsFile, CASEINSEN, "F_SAMP", &colnum, &status);
  sprintf( errorMessage, "Could not find f_samp column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  fits_read_col( thisFitsFile, TDOUBLE, colnum, (foundDet+1), 1, 1, NULL, fsample, NULL, &status);
  sprintf( errorMessage, "Could not read f_samp for detector %s in file %s", formatStruct.detectorID, theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  fits_get_colnum( thisFitsFile, CASEINSEN, "NET_RJ", &colnum, &status);
  sprintf( errorMessage, "Could not find net_rj column in file %s", theFile->name );
  M3_ErrorCheck( -1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  fits_read_col( thisFitsFile, TDOUBLE, colnum, (foundDet+1), 1, 1, NULL, NET, NULL, &status);
  sprintf( errorMessage, "Could not read net_rj for detector %s in file %s", formatStruct.detectorID, theFile->name );
  M3_ErrorCheck( -1, errorMessage, detids, M3_EFLAG_DEFAULT );

  /* close the file */

  fits_close_file( thisFitsFile, &status);
  sprintf( errorMessage, "Error closing fits file named %s", theFile->name );
  M3_ErrorCheck(-1, errorMessage, status == 0, M3_EFLAG_DEFAULT );

  return;
}


void M3_CreateNoiseFilter( double fknee, double alpha, double fmin, double fsample, double NET, int64_t corr, double *itt_noise )
{
  /* parameters */
  char errorMessage[M3_CONTENT_LENGTH];
  char errorString[M3_CONTENT_LENGTH];
  char strBuffer[M3_CONTENT_LENGTH];
  char strBuffer2[M3_CONTENT_LENGTH];
  char invnttfile[M3_CONTENT_LENGTH];
  char invntttext[M3_CONTENT_LENGTH];
  char psdfile[M3_CONTENT_LENGTH];
  int apo = 0;
  int apowidth = 0;
  double rttwo = sqrt(2.0);

  /* variables */
  int i;  
  int tfields;
  int npsd;
  int nfft;
  double *freq;
  double *psd;
  double *fft;
  double *invntt;
  
  double nyquist;
  double bandwidth;
  double halfband;
  double cur;
  fftw_plan plan;

  /* temp variables */

  double temp;
  double ktemp;
  double mtemp;
  char command[M3_CONTENT_LENGTH];


  if (corr == 0) {
    itt_noise[0] = 1.0 / (NET * NET * fsample);
    return;
  }

  nfft = 2;
  while (nfft < 2*corr) {
    nfft *= 2;
  }
  nfft *= M3_NOISE_GEN_CORRFACTOR;
  
  nfft = (nfft < M3_NOISE_GEN_MINIMUM) ? M3_NOISE_GEN_MINIMUM : nfft;

  /* generate the psd */

  npsd = 1 + (int)(nfft / 2);
  nyquist = (int)(fsample / 2);
  bandwidth = fsample / (double)(nfft);
  halfband = 0.5 * bandwidth;
  
  freq = (double*)calloc(npsd, sizeof(double));
  psd = (double*)calloc(npsd, sizeof(double));
  M3_ErrorCheck(-1, "M3_File_readData_LevelS_noise()", (freq && psd), M3_EFLAG_MALLOC_FSTRING );

  ktemp = pow(fknee, alpha);
  mtemp = pow(fmin, alpha);  
  freq[0] = 0.0;
  psd[0] = NET * sqrt( ktemp / mtemp );
  for (i = 1; i < npsd-1; i++) {
    cur = (double)i * bandwidth;
    freq[i] = cur;
    temp = pow(cur, alpha);
    psd[i] = rttwo * NET * sqrt((temp + ktemp) / (temp + mtemp));
  }
  freq[npsd-1] = bandwidth * (double)(npsd-1);
  temp = pow(freq[npsd-1], alpha);
  psd[npsd-1] = NET * sqrt((temp + ktemp) / (temp + mtemp));
  
  /* convert PSD to (N / |Y(f)|^2) and pre-normalize for FFTW */

  invntt = (double*)calloc(nfft, sizeof(double));
  fft = (double*)calloc(nfft, sizeof(double));
  M3_ErrorCheck(-1, "M3_File_readData_LevelS_noise()", (invntt && fft), M3_EFLAG_MALLOC_FSTRING );
  
  plan = fftw_plan_r2r_1d(nfft, fft, invntt, FFTW_HC2R, FFTW_ESTIMATE);
  
  temp = (double)(nfft) * (double)(nfft);
  for (i = 0; i < npsd; i++) {
    if ((i > 0) && (i < npsd-1)) {
      fft[i] = 1.0 / (temp * halfband * (psd[i] * psd[i]));
    } else {
      fft[i] = 1.0 / (temp * bandwidth * (psd[i] * psd[i]));
    }
  }
  
  /* compute inverse noise correlations.  */
  fftw_execute(plan);
  memcpy( itt_noise, invntt, (corr+1)*sizeof(double));

  fftw_destroy_plan(plan);
  fftw_cleanup();
  free(psd);
  free(invntt);
  free(fft);
  free(freq);

  return;
}


void M3_File_readData_LevelS_noise( M3_File *theFile,  double *itt_noise )
{
  double fsample = 0.0;
  double fknee = 0.0;
  double fmin = 0.0;
  double NET = 0.0;
  double alpha = 0.0;
  int64_t corr = 0.0;

  M3_File_getFPDBparam( theFile, &fknee, &alpha, &fmin, &fsample, &NET );
  corr = theFile->param.noise.corLength;
  M3_CreateNoiseFilter( fknee, alpha, fmin, fsample, NET, corr, itt_noise );
}

---