/**************************************************************************/
/* makemap.c - generate a HEALPix map with hpic                           */
/**************************************************************************/

#include <stdio.h>
#include "math.h"
#include "hpic.h"

float spherefunction(double theta, double phi);

int main(int argc, char **argv)
{  
  char *outputfilename;

  /* check input */
  if (argc!=2) {
    printf("usage: makemap myoutputmap.fits\n");
    exit(EXIT_FAILURE);
  }
  outputfilename = argv[1];  
  if (outputfilename == NULL) {
    printf("usage: makemap myoutputmap.fits\n");
    exit(EXIT_FAILURE);
  }

  size_t nmaps = 1;
  size_t nside = 128; //this MUST be a power of 2!
  int order = HPIC_NEST;
  int coord = HPIC_COORD_C;

  char filename[HPIC_STRNL];
  char extname[HPIC_STRNL];
  char comment[HPIC_STRNL];
  char creator[HPIC_STRNL];
  strncpy(filename,outputfilename,HPIC_STRNL);
  strncpy(extname,"extname",HPIC_STRNL);
  strncpy(creator,"makemap",HPIC_STRNL);
  strncpy(comment,"my test map",HPIC_STRNL);

  hpic_fltarr* hpic_maps; 			
  hpic_float* hpic_map;
  hpic_maps = hpic_fltarr_alloc(nmaps);
  hpic_map = hpic_float_alloc(nside,order,coord,HPIC_STND);
  hpic_fltarr_set(hpic_maps,0,hpic_map);

  size_t i;
  size_t npix = 12*nside*nside;
  float value;
  double theta, phi;
  
  for(i=0;i<npix;i++){
    hpic_pix2ang_nest(nside,i,&theta,&phi);
    value = spherefunction(theta,phi);
    hpic_float_set(hpic_map,i,value);
  }

  hpic_keys *keys = hpic_keys_alloc();
  hpic_fits_full_write(filename,creator,extname,comment,hpic_maps,keys);

  return EXIT_SUCCESS;
}

/* any old function on the sphere (theta, phi in radians) */
float spherefunction(double theta, double phi)
{  
  return (float)cos(8.0*theta)*sin(5.0*phi);
}

/**************************************************************************/