M3: binMap.c

/*******************************************************************************
*   M3:  binMap.c                                                              *
*                                                                              *
*   Version 2.0.2 May 2007                                                     *
*                                                                              *
*   Copyright (C) 2004-2007  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  *
*                                                                              *
*******************************************************************************/


/*  Christopher M. Cantalupo November 2007 */

#include <stdlib.h>
#include <M3.h>

#define STRING_LENGTH 256

int main( int argc, char **argv )
{
  M3_RunConfigStruct *runConfig;
  M3_DataSetLL *dataSet;
  M3_IntervalLL *intervalNode;
  int notDone;
  M3_PointingEl *pointing = NULL;
  M3_PixelClassLL *pixelClass;
  M3_PointingClassLL *pointingClass;
  int32_t numNZ, numClassNZ, numPixelClass, mapSize;
  M3_MapEl **mapArray;
  double **weightArray;
  int32_t *nnzArray;
  int32_t numPixel;
  int64_t i, j, k, l, n;
  M3_Interval interval;
  double *tod = NULL;
  char mapName[STRING_LENGTH];
  FILE *fid;
  M3_MapHeader mapHeader;
  char temp = 0;

  /* Read run config from disk */
  M3_RunConfigStruct_Read(&runConfig, argv[1], M3_SPRINGTIDE_RUN_TYPE );

  /* Count the number of pixel classes */
  numPixelClass = 0;
  for( notDone = M3_RunConfigStruct_GetPixelClassRoot(runConfig, &pixelClass); notDone; notDone = M3_PixelClassLL_GetNextNodeInList(&pixelClass) )
    numPixelClass++;

  /* Allocate some space */
  mapArray = (M3_MapEl **)calloc(numPixelClass, sizeof(M3_MapEl *));
  weightArray = (double **)calloc(numPixelClass*numPixelClass, sizeof(double *));
  nnzArray = (int32_t *)malloc(numPixelClass*sizeof(int32_t));

  /* Allocate maps and weights for each pixel class */
  for( i = 0, notDone = M3_RunConfigStruct_GetPixelClassRoot(runConfig, &pixelClass); 
       notDone; 
       i++, notDone = M3_PixelClassLL_GetNextNodeInList(&pixelClass) )
  {
    M3_PixelClassLL_GetNumPixelInClass( pixelClass, &mapSize);
    mapArray[i] = (M3_MapEl *)calloc(mapSize,sizeof(M3_MapEl));
    /* Initialize map */
    for( j = 0; j < mapSize; j++ )
    {  
      mapArray[i][j].pixel = j;
      mapArray[i][j].value = nan(&temp);
    }
    weightArray[i] = calloc(mapSize,sizeof(double));
  }


  /* Go through all of the dataSets (detectors) */
  for ( notDone = M3_RunConfigStruct_GetDataSetRoot( runConfig, &dataSet );
        notDone;
        M3_DataSetLL_GetNextNodeInList( &dataSet ))
  {
    /* Figure out how many non-zeros per row of the pointing matrix */
    M3_DataSetLL_GetNumNZ(dataSet, &numNZ);

    /* Go through the covered intervals */
    for( notDone = M3_DataSetLL_GetCoveredIntervalRoot( dataSet, &intervalNode);
         notDone;
         notDone = M3_IntervalLL_GetNextNodeInList( &intervalNode ) )
    {
      /* Retrieve the interval structure from the list node */
      M3_IntervalLL_GetInterval( intervalNode, &interval );
      /* Count the number of samples in the interval*/
      n = interval.lastSample - interval.firstSample + 1;

      /* If we have alread allocated tod and pointing, call realloc */
      if( tod && pointing )
      {
        tod = (double*)realloc( tod, sizeof(double)*n);
        pointing = (M3_PointingEl*)realloc( pointing, numNZ*sizeof(M3_PointingEl)*n );
      }
      /* Otherwise allocate memory with malloc */
      else
      {
        tod = (double*)malloc( sizeof(double)*n);
        pointing = (M3_PointingEl*)malloc( numNZ*sizeof(M3_PointingEl)*n );
      }

      /* Get the pointing data and the time ordered data*/
      M3_DataSetLL_GetPointing(dataSet, interval, pointing );
      M3_DataSetLL_GetTOD( dataSet, interval, tod );

      /* Cycle over all of the samples in the interval.   
         l keeps track of where you are in the pointing array */
      l = 0;
      /* i keeps track of the time index from the beginning of the interval */
      for( i = 0; i < n; i++ )
      {
        /* j keeps track of the pixel class index */
        j = 0;
        /* Get the root pointing class for this data set */
        notDone = M3_DataSetLL_GetPointingClassRoot( dataSet, &pointingClass );
        /* Go over the poiting classes */
        while( notDone )
        {
          /* Find out how many non-zeros per time sample this poining class contributes */
          M3_PointingClassLL_GetNumNZ( pointingClass, &numClassNZ );
          /* Loop over the non-zeros */
          for( k = 0; k < numClassNZ; k++ )
          {
            /* If the map element has already been set once, we can simply add to the current value */
            if(!isnan(mapArray[j][pointing[l].pixel].value))
            {
              /* Sum the weighted TOD into map */
              mapArray[j][pointing[l].pixel].value += tod[i]*pointing[l].weight;
              /* Sum the square of the weights */
              weightArray[j][pointing[l].pixel] += pointing[l].weight*pointing[l].weight;
            }
            /* Otherwise we need to set the value.  */
            else
            {
              /* Set the map to the weighted TOD value */
              mapArray[j][pointing[l].pixel].value = tod[i]*pointing[l].weight;
              /* Set the weight to the square of the weights */
              weightArray[j][pointing[l].pixel] = pointing[l].weight*pointing[l].weight;
            }
            l++;
          }
          j++;
          /* Go to the next pointing class */
          notDone = M3_PointingClassLL_GetNextNodeInList( &pointingClass ); 
        }
      }
    }

    /* Go over the pixel classes */
    M3_RunConfigStruct_GetPixelClassRoot( runConfig, &pixelClass );
    for( i = 0; i < numPixelClass; i++ )
    {
      /* Now divide the map by the sum of the squares of the weights 
         note that we have not actually solved for Q and U, but are 
         mearly finding the average weighted value.  */
      M3_PixelClassLL_GetNumPixelInClass( pixelClass, &numPixel );
      k = 0;
      for( j = 0 ; j < numPixel; j++)
      {
        if( !isnan(mapArray[i][j].value) && weightArray[i][j] )
        {
          mapArray[i][k].value = mapArray[i][j].value/sqrt(weightArray[i][j]);
          mapArray[i][k].pixel = mapArray[i][j].pixel;
          k++;
        }
      }
      /* Get the name of the map for output */
      M3_PixelClassLL_GetMapName( pixelClass, mapName, STRING_LENGTH );
      /* Write the map to disk in MADCAP3 format */
      fid = fopen( mapName, "w");
      mapHeader.numPixelInClass = numPixel;
      mapHeader.numPixelInMap = k;
      fwrite( &mapHeader, 1, sizeof(M3_MapHeader), fid);
      fwrite( mapArray[i], sizeof(M3_MapEl), k, fid);
      fclose(fid);
      /* Go to the next pixel class */
      M3_PixelClassLL_GetNextNodeInList( &pixelClass );
    }
  }

  free(tod);
  free(pointing);
  return 0;
}