#include <stdlib.h>
#include <math.h>
#include "ph_dcm.h"

#define MAX_30_BINS 5
#define MAX_PH_BINS 121
#define MAX_PH_SAMPLES 1890

#if 1
static double theta[24] = {    /*   -90 < theta < 90   latitude */
  -78.7500, -56.2500, -39.3750, -28.1250,
  -19.6875, -14.0625, -8.43750, -2.81250,
   2.81250,  8.43750,  14.0625,  19.6875,
   28.1250,  39.3750,  56.2500,  78.7500,
   78.7500,  56.2500,  33.7500,  11.2500, /* back half */
  -11.2500, -33.7500, -56.2500, -78.7500};
#else
static double theta[24] = {    /*   0 < theta < 180   Co-latitude */
  168.75000, 146.25000, 129.37500, 118.12500,
  109.68750, 104.06250, 98.437500, 92.812500,
  87.187500, 81.562500, 75.937500, 70.312500,
  61.875000, 50.625000, 33.750000, 11.250000,
  168.75000, 146.25000, 123.75000, 101.25000, /* back half */
  78.750000, 56.250000, 33.750000, 11.250000};
#endif

static double dtheta[24] = {22.5,22.5,11.25,11.25,5.625,5.625,5.625,5.625,
			    5.625,5.625,5.625,5.625,11.25,11.25,22.5,22.5,
			    22.5,22.5,22.5,22.5,22.5,22.5,22.5,22.5};
static double phi[32] = {
   174.37500,  163.12500,  151.87500,  140.62500,  129.37500,  118.12500,
   106.87500,  95.625000,  84.375000,  73.125000,  61.875000,  50.625000,
   39.375000,  28.125000,  16.875000,  5.6250000, -5.6250000, -16.875000,
  -28.125000, -39.375000, -50.625000, -61.875000, -73.125000, -84.375000,
  -95.625000, -106.87500, -118.12500, -129.37500, -140.62500, -151.87500,
  -163.12500, -174.37500 };
static double dphi = 11.25;
static double domegas[24] = {
  0.014946219,  0.042563230,  0.029754132,  0.033946244,
  0.018142453,  0.018691368,  0.019060275,  0.019245620,
  0.019245620,  0.019060275,  0.018691368,  0.018142453,
  0.033946244,  0.029754132,  0.042563230,  0.014946219,
  0.014946219,  0.042563230,  0.063700375,  0.075139715, /* back half.  added by fvm */
  0.075139715,  0.063700375,  0.042563230,  0.014946219};
static float geom_area_ph[24] = {4,4,2,2,1,1,1,1,1,1,1,1,2,2,4,4,4,4,4,4,4,4,4,4};
static int tsect[32] = {0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6, 7,
                        8, 9,10,11,12,12,13,13,14,14,14,14,15,15,15,15};

static int no_blank[32]   = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
			     0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
static int all_blank[32]  = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
			     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
static int back_blank[24] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
			     1,1,1,1,1,1,1,1}; /*back half of detector blanked*/


int get_ph_mapcode_idl(int argc,void *argv[])
{
  idl_ph_pk_data idl;
  static data_map_3d map;                 
  packet_selector *pksp;
  static packet_selector pks_ph;
  uint4 id;
  static uint2 spin;                      
  packet *pk;

  idl.options = (int4   *)argv[0];  /* get input args */
  idl.time    = (double *)argv[1];
  idl.advance = (int2   *)argv[2];
  idl.mapcode = (int4   *)argv[3];
  idl.idtype  = (uint2  *)argv[4];
  idl.instseq = (uint2  *)argv[5];

  pksp = &pks_ph;    
  id = P3D_ID;

  /* negitive option[1] means get by time*/
  if     (*idl.advance)    {SET_PKS_BY_INDEX(*pksp,pksp->index+*idl.advance,id);}
  else if(idl.options[1]<0){SET_PKS_BY_TIME (*pksp,idl.time[0],id);             }
  else                     {SET_PKS_BY_INDEX(*pksp,idl.options[1],id);          }

  pk = get_packet(pksp);             
  if (pk==NULL) return(0);

  *idl.time    = pk->time;
  *idl.idtype  = pk->idtype;
  *idl.instseq = pk->instseq;

  idl.mapcode[0] = get_ph_map_code(pk);
 
  return(1);
}

int get_phb_mapcode_idl(int argc,void *argv[])
{
  idl_ph_pk_data idl;
  static data_map_3d map;                 
  packet_selector *pksp;
  static packet_selector pks_ph;
  uint4 id;
  static uint2 spin;                      
  packet *pk;

  idl.options = (int4   *)argv[0];  /* get input args */
  idl.time    = (double *)argv[1];
  idl.advance = (int2   *)argv[2];
  idl.mapcode = (int4   *)argv[3];
  idl.idtype  = (uint2  *)argv[4];
  idl.instseq = (uint2  *)argv[5];

  pksp = &pks_ph;    
  id = P3D_BRST_ID;

  /* negitive option[1] means get by time*/
  if     (*idl.advance)    {SET_PKS_BY_INDEX(*pksp,pksp->index+*idl.advance,id);}
  else if(idl.options[1]<0){SET_PKS_BY_TIME (*pksp,idl.time[0],id);             }
  else                     {SET_PKS_BY_INDEX(*pksp,idl.options[1],id);          }

  pk = get_packet(pksp);             
  if (pk==NULL) return(0);

  *idl.time    = pk->time;
  *idl.idtype  = pk->idtype;
  *idl.instseq = pk->instseq;

  idl.mapcode[0] = get_ph_map_code(pk); 
  return(1);
}

int get_next_ph5_struct(packet_selector *pks, idl_ph5 *dat5){

  static uint2 spin;
  packet *pk;
  pk = get_packet(pks);
  
  if(pk) spin = pk->spin;
  
  /* gather all packets of this type in this spin */
  while(pk){
    if(pk->quality & (~pkquality)){
	dat5->time = pk->time;
	return(0);
    }
    decom_map_ph5(pk,dat5);
    spin = pk->spin;
    if (!pk->next || spin != (pk->next)->spin) break ;
    pk = pk->next;
  }
  return(pk ? 1 : 0);  
}
    
int get_next_ph97_struct(packet_selector *pks, idl_ph97 *dat97)
{
  static uint2 spin;
  packet *pk;
  pk = get_packet(pks);
  
  if(pk) spin = pk->spin;
  /* gather all packets of this type in this spin */

  while(pk){
    if(pk->quality & (~pkquality)){
	dat97->time = pk->time;
	return(0);
    }
    decom_map_ph97(pk,dat97);
    spin = pk->spin;
    if (!pk->next || spin != (pk->next)->spin) break ;
    pk = pk->next;
  }
  return(pk ? 1 : 0);
}

int get_next_ph121_struct(packet_selector *pks, idl_ph121 *dat121){

  static uint2 spin;
  packet *pk;
  pk = get_packet(pks);
 
  if(pk) spin = pk->spin;
  
  /* gather all packets of this type in this spin */
  while(pk){
    if(pk->quality & (~pkquality)){
	dat121->time = pk->time;
	dat121->valid = 0;
	return(0);
    }
    decom_map_ph121(pk,dat121);
    spin = pk->spin;
    if (!pk->next || spin != (pk->next)->spin) break ;
    pk = pk->next;
  }
  return(pk ? 1 : 0);  
}

int get_next_ph56_struct(packet_selector *pks, idl_ph56 *dat56)
{
  static uint2 spin;
  packet *pk;
  pk = get_packet(pks);
  
  if(pk) spin = pk->spin;
  /* gather all packets of this type in this spin */

  while(pk){
    if(pk->quality & (~pkquality)){
	dat56->time = pk->time;
	return(0);
    }
    decom_map_ph56(pk,dat56);
    spin = pk->spin;
    if (!pk->next || spin != (pk->next)->spin) break ;
    pk = pk->next;
  }
  return(pk ? 1 : 0);
}

int get_next_ph65_struct(packet_selector *pks, idl_ph65 *dat65)
{
  static uint2 spin;
  packet *pk;
  pk = get_packet(pks);
  
  if(pk) spin = pk->spin;
  /* gather all packets of this type in this spin */

  while(pk){
    if(pk->quality & (~pkquality)){
	dat65->time = pk->time;
	return(0);
    }
    decom_map_ph65(pk,dat65);
    spin = pk->spin;
    if (!pk->next || spin != (pk->next)->spin) break ;
    pk = pk->next;
  }
  return(pk ? 1 : 0);
}

int get_next_ph88_struct(packet_selector *pks, idl_ph88 *dat88)
{
  static uint2 spin;
  packet *pk;
  pk = get_packet(pks);
  
  if(pk) spin = pk->spin;
  /* gather all packets of this type in this spin */

  while(pk){
    if(pk->quality & (~pkquality)){
	dat88->time = pk->time;
	return(0);
    }
    decom_map_ph88(pk,dat88);
    spin = pk->spin;
    if (!pk->next || spin != (pk->next)->spin) break ;
    pk = pk->next;
  }
  return(pk ? 1 : 0);
}

/* This subroutine will sort an uint array of length arraylength
   It is not the fastest way, but with typically 5 elements to sort,
   it is fairly quick */
void bubble_sort(uint array[], int arraylength){ /*modified bubble sort*/
  int i,j,swap;
  uint temp;
  for(i=arraylength-1;i>=0;--i){
    swap=0;
    for(j=0;j<i;j++){
      if(array[j]>array[j+1]){
	temp       = array[j];
	array[j]   = array[j+1];
	array[j+1] = temp;
	swap=1;
      }
    }
    if(!swap) return;
  }
  return;
}

/* This subroutine searches the ptmap for the bins in ptmap
   which have 30 energy steps. */
int find_ph_e30bins(uint ptmap[32*32],uint e30bins[MAX_30_BINS]){
  int i,e30,e30b,got_it,nfound;
  i=e30=e30b=got_it=nfound=0;

  while(nfound<MAX_30_BINS){
    /* find next e30 bin, or quit when entire array searched */
    while((ptmap[i] & 0xc000) != E30) { if(++i==24*32) break; }  
    e30b=(ptmap[i++] & 0x3fff);                 /*get the bin number*/
    if (nfound==0) {                            /*if first e30 bin*/
      e30bins[nfound++]=e30b;                     /*save the bin number*/
    } else {                                    /*else*/
      got_it=0;
      for(e30=0;e30<nfound;e30++)                 /*see if bin found already*/
	if (e30bins[e30]==e30b) got_it++;
      if(!got_it) e30bins[nfound++]=e30b;         /*if not found, save it */
    }
    if (nfound==MAX_30_BINS) break;
  }
  bubble_sort(e30bins,MAX_30_BINS);           /*sort the bin array*/
  return(nfound);
}

/* this subroutine fills the 30 energy step bin data structures
   for a (non-burst) ph map. */
int decom_map_ph5(packet *pk, idl_ph5 *map){
  int i,offset,seq;
  static int status,nsamples,npackets;
  int ndata,max_bytes_per_packet=390;
  uchar shift,*data;
  int b,b30,e,p,t;
  static BinRange bin[MAX_30_BINS]; 
  static float alldata[MAX_30_BINS*30];
  double dt,nrg[30],dnrg[30];
  static uint *ptmap, e30bins[MAX_30_BINS];
  uchar p0=28;
  static int num_phi[MAX_30_BINS];  /* moved up from next block */

  i=offset=seq=b=b30=e=p=t=shift=ndata=npackets=0;

  if(pk->spin != map->spin || status<=0){
    /* This is executed only on the first of "npackets" calls.
       It gets all of the angular and energy information, and 
       sets things such as the time and spin. */
    PCFG *pcfg;   
    static int theta_slots[32][MAX_30_BINS];
    static int phi_slots[24/*16*/][MAX_30_BINS];
    static int num_theta[MAX_30_BINS];

    /* set to NaN instead of 0 because the last bin may be turned off
       this way the data will be ignored instead of falsly represented as 0 */
    for (i=0;i<MAX_30_BINS*30;i++) alldata[i] = NaN;

    switch(map->mapcode){
    case(MAP11d):  nsamples=1890;  npackets = 5;  ptmap=ptmap_5;    break;
    case(MAP11b):  nsamples=1530;  npackets = 4;  ptmap=ptmap_11b;  break;
    default:
      printf("decom_map_ph5 (in ph_dcm.c): unknown mapcode: %d\n\r",
	     map->mapcode);
      return(0);
    }
    find_ph_e30bins(ptmap,e30bins); /*find the bins with 30 energy steps*/
                                    /*eg: for MAP11d: e30bins={0,16,35,51,120}*/
    pcfg = get_PCFG(pk->time);

    map->time        = pk->time;
    map->integ_t     = get_3d_integ_t(pk) * pcfg->spin_period;
    map->delta_t     = get_3d_delta_t(pk) * pcfg->spin_period;
    map->end_time    = map->time + map->integ_t;
    map->geom_factor = pcfg->ph_geom;
    map->spin        = pk->spin;
    /* these are set in idl
    map->mass        = 1836*5.6856591e-6;
    map->nbins       = MAX_30_BINS;
    map->nenergy     = 30; */

    dt       = map->integ_t / 32 / 32; /* nspins*spin period/nphis/nenergies */
    shift    = (pk->instseq >> 3) & 0x1f;
    map->shift = shift;
    status   = npackets;  /*decrement variable, counts packets left to gather*/

    for(i=0;i<30*4;i++){
      map->dac_code[i] = pcfg->phdac_tbl[i];
      map->volts[i]    = pcfg->ph_volts_tbl[i];
    }
    
    for(b=0;b<map->nbins;b++){  /* zero the static struct "bin" */
      bin[b].ne      = 0;       /* and a few other binning aids */
      bin[b].offset  = 0;
      bin[b].geom    = 0.;
      bin[b].theta   = 0.;
      bin[b].dtheta  = 0.;
      bin[b].dphi    = 0.;
      bin[b].domega  = 0;
      num_theta  [b] = 0;
      num_phi    [b] = 0;
      map->domega[b] = 0.;
      for(e=0;e<map->nenergy;e++){
        map->dtheta[b][e] = 0.;
        map->dphi  [b][e] = 0.;
	map->nrg  [b][e] = 0.;
	map->dnrg [b][e] = 0.;
	map->theta[b][e] = 0.;
	map->phi  [b][e] = 0.;
	map->flux [b][e] = 0.;
	map->gf   [b][e] = 0.;
	map->dt   [b][e] = 0.;
      }
      for (e=0; e<30; e++) bin[b].phi[e]     = 0.;
      for (p=0; p<32; p++) theta_slots[p][b] = 0;
      for (t=0; t<24; t++) phi_slots[t][b]   = 0;
    }
    
    for(p=0;p<32;p++){
      for(t=0;t<24;t++){
	
	b = ptmap[t*32+p];        /* get the bin number for a phi/theta pair */
	if ((b & 0xc000) != E30) continue;     /*we only want high res. bins */
	b &= 0x3fff;                       /* chop off the energy steps info */
	
	b30=0;
	while (b30<MAX_30_BINS)        /* find b30 such that b==e30bins[b30] */
	  if (b==e30bins[b30]) break; else b30++;

	if(b30 > map->nbins){                     
	  fprintf(stderr,"3D bin error\n\r");
	  return(0);
	}

	bin[b30].ne=30;

	/*if(t<16){*/
	bin[b30].phi[0]     += phi[p];
	bin[b30].theta      += theta[t] * dtheta[t];
	bin[b30].dphi       += dphi;
	bin[b30].dtheta     += dtheta[t];
	bin[b30].domega     += domegas[t];
	bin[b30].offset     += 1;
	theta_slots[p][b30] += 1;
	phi_slots[t][b30]   += 1;
	/*}*/
	bin[b30].geom       += geom_area_ph[t];
      }
    }
    for(b=0;b<map->nbins;b++){
      float phi0 = bin[b].phi[0];
      for (p=0; p<32; p++)
	if (num_theta[b] < theta_slots[p][b])
	  num_theta[b] = theta_slots[p][b];
      for (t=0; t<24/*16*/; t++)
	if (num_phi[b] < phi_slots[t][b])
	  num_phi[b] = phi_slots[t][b];
      if(bin[b].offset ==0){
	/*fprintf(stderr,"Invalid bin #%d\n\r",b); */
	continue;
      }
      /*printf("bin[%1d].offset: %4d\n\r",b,bin[b].offset);*/
      bin[b].theta /= bin[b].dtheta;
      phi0   /= bin[b].offset;
      phi0   += ((float)shift + (float)p0) * 360./32.;
      if(phi0 < 0.)
	phi0 +=360.;
      if(phi0 >= 360.)
	phi0 -=360.;
      bin[b].dtheta /= num_phi[b];
      bin[b].dphi   /= num_theta[b];
      bin[b].geom   /= num_phi[b];    
      
      /* build up phi across energy bins */
      for (e=0; e<bin[b].ne; e++)
	bin[b].phi[e] = phi0 +	(.5 - (e+.5)/bin[b].ne)*dphi;
    }
 
/* this may need fixing.  the question is do i make it compatible with
   maps that have ne=0 bins?  that would add another subroutine or make
   this block considerably larger.  currently, no ph maps have bins with ne=0*/
    /* the following solution is not compatible with ne=0 bins in ptmap*/
    offset=0;
    for(b=0;b<map->nbins;b++){
      if (b==0) bin[b].offset=e30bins[0]*15;
      else      bin[b].offset=bin[b-1].offset+(e30bins[b]-e30bins[b-1]-1)*15+30;
    }

    get_esteps_ph(nrg,map->nenergy,MIDDLE,pcfg);
    get_esteps_ph(dnrg,map->nenergy,WIDTH,pcfg);

  }

  seq    = pk->instseq & 0x07;
  ndata  = max_bytes_per_packet;
  offset = seq * ndata;
  if (ndata>(int)pk->dsize) ndata = pk->dsize;
  data   = pk->data;


  for(b=0;b<map->nbins;b++){
    if(bin[b].offset <  offset)          continue;  /*data is in earlier pk*/
    if(bin[b].offset >  ndata+offset-30) break;     /*data is in later   pk*/
    if(bin[b].offset >  nsamples-30){               /*offset too large     */
      char instStr[12]="PESA HIGH";
      fprintf(stderr,"3D binning error: offset(%d) >= nsamples(%d)\n\r",
	      bin[b].offset+30,nsamples);
      fprintf(stderr,"    ndata: %d, map type: %s\n\r", ndata, instStr);
    } else {
      for(e=0;e<30;e++) 
	alldata[b*30+e]=decomp19_8(*(data+bin[b].offset-offset+e));
    }
  }
  status--;
  
  if(status == 0){   /* if all packets received, fill the data structure */
    for(b=0;b<map->nbins;b++){
      map->domega[b] = bin[b].domega;
      for(e=0;e<map->nenergy;e++){
        map->dtheta[b][e] = bin[b].dtheta;
        map->dphi  [b][e] = bin[b].dphi;
	map->nrg  [b][e] = nrg[e];
	map->dnrg [b][e] = dnrg[e];
	map->theta[b][e] = bin[b].theta;
	map->phi  [b][e] = bin[b].phi[e];
	map->gf   [b][e] = bin[b].geom;
	map->dt   [b][e] = num_phi[b] * dt;
	map->flux [b][e] = alldata[b*30+e];
	/*printf("bin: %d: e: %d, ne: %d, offset: %d, data: %f\n\r",
	       b,e,bin[b].ne,b*30+e,alldata[b*30+e]);*/
      }
    }

    map->valid=1;
  }
  return(1);  
}

/* This subroutine gets the angular information for the low res data structures
   for decom_map_ph{97,121} (others may be added with minimal effort) */
int make_ph_ptmap(PCFG *pcfg,uchar shift,int4 mapcode,int num_phi[MAX_PH_BINS],
		 uint ptmap[32*32],BinRange bin[]){
  int b,e,offset,ne,p,t,nsamples,nbins;
  int *t_blank, *p_blank;
  static int first_call=1;
  uint *local_ptmap,ps=0;
  uchar p0=28;
  static int theta_slots[32][MAX_PH_BINS],phi_slots[16][MAX_PH_BINS];
  static int num_theta[MAX_PH_BINS];
  /*static int num_phi[MAX_PH_BINS]; made into function argument,
    used by calling functions*/
  
  static struct {     /* if any of these value change, must recalculate map */
    uchar   shift;                      /* shift between p and phi in ptmap */
    int4    mapcode;                                            /* data map */
    nvector esa_swp_select;                        /* ummm, not really sure */
    uint2   select_sector;                                         /* ditto */
    uchar   esa_pha_basech;                                        /* ditto */
    nvector bts_c_val;                                             /* ditto */
  } prev_mapinfo, mapinfo;            /* intial values are bogus */

  mapinfo.shift          = shift;
  mapinfo.mapcode        = mapcode;
  mapinfo.esa_swp_select = pcfg->norm_cfg.esa_swp_select;
  mapinfo.select_sector  = pcfg->norm_cfg.select_sector; 
  mapinfo.esa_pha_basech = pcfg->norm_cfg.esa_pha_basech;
  mapinfo.bts_c_val      = pcfg->norm_cfg.bts_c_val;
	    
/*What are the nec. and suf. conditions such that all of 
  the angular info will be the same this time as last call?
  0: make_ph_ptmap                   must not be first call to make_ph_ptmap
  1: mapcode                         must be the same
  2: shift                           must be the same 
  4: pk->time and last_ph_time are on different sides of PH_BACK_OFF_TIME */
	
  if (!first_call) 
    if (memcmp(&prev_mapinfo,&mapinfo,sizeof(mapinfo))==0)
      return(1);  /* packet info is the same. use old angular info */
  
  first_call=0;

  switch (mapinfo.bts_c_val >> 14) {
  case 3:
    t_blank = no_blank;
    p_blank = no_blank;
    break;
  case 2:
    t_blank = back_blank;  /* blank all phis/theta pairs */
    p_blank = all_blank;   /* with back facing thetas */
    break;
  default:
    printf("make_ph_ptmap: bts_c_val unknown: %x\r\n",mapinfo.bts_c_val);
    return(0);
  }
	

  switch(mapcode){
  case(MAP11d): nbins=121;  nsamples=1890;  local_ptmap=ptmap_5;    break;
  case(MAP11b): nbins= 97;  nsamples=1530;  local_ptmap=ptmap_11b;  break;
  case(MAP_8):  nbins= 56;  nsamples= 960;  local_ptmap=ptmap_8;    break;
  case(MAP_0):  nbins= 65;  nsamples=1050;  local_ptmap=ptmap_0;    break;
  case(MAP22d): nbins= 88;  nsamples=1380;  local_ptmap=ptmap_22d;   break;
  default:
    printf("make_ph_ptmap: map index unknown: %d\n",mapcode);
    return(0);
  }
  
  for(b=0;b<nbins;b++){  /* zero the static struct "bin" */
    bin[b].ne      = 0;       /* and a few other binning aids */
    bin[b].offset  = 0;
    bin[b].geom    = 0.;
    bin[b].theta   = 0.;
    bin[b].dtheta  = 0.;
    bin[b].dphi    = 0.;
    bin[b].domega  = 0;
    num_theta  [b] = 0;
    num_phi    [b] = 0;
    for (e=0; e<30; e++) bin[b].phi[e]     = 0.;
    for (p=0; p<32; p++) theta_slots[p][b] = 0;
    for (t=0; t<16; t++) phi_slots  [t][b] = 0;
  }
  
  for(p=0;p<32;p++){
    for(t=0;t<24;t++){
      b = local_ptmap[t*32+p]; /*get the bin number for a phi/theta pair*/

      switch (b & 0xc000) {    /*get the # of energy steps for this bin */ 
      case E30:  ne = 30;  break;
      case E0:   ne =  0;  break;
      default:   ne = 15;  break;
      }
      
      b &= 0x3fff;           /* chop off the energy steps info */
      
      if(b > nbins && ne ){  /* check order of operations */
	fprintf(stderr,"3D bin error\n\r");
	return(0);
      }
      
      ps = (p - shift) & 0x1f;
      ptmap[t*32+ps] = (ne == 0) ? -1 : b;
      if(ne == 0)
	continue;
      if (b < nbins ) bin[b].ne = ne;
      
      if(t<16){
	bin[b].phi[0]     += phi[p];
	bin[b].theta      += theta[t] * dtheta[t];
	bin[b].dphi       += dphi;
	bin[b].dtheta     += dtheta[t];
	bin[b].domega     += domegas[t];
	bin[b].offset     += 1;
	theta_slots[p][b] += 1;
	phi_slots[t][b]   += 1;
      }
      if(p_blank[p] & t_blank[t])
	continue;     /* skip blanked anodes */
      bin[b].geom += geom_area_ph[t];
    }
  }

  for(b=0;b<nbins;b++){
    float phi0 = bin[b].phi[0];
    for (p=0; p<32; p++)
      if (num_theta[b] < theta_slots[p][b]) num_theta[b] = theta_slots[p][b];
    for (t=0; t<16; t++)
      if (num_phi[b]   < phi_slots[t][b])   num_phi[b]   = phi_slots[t][b];
    if(bin[b].offset==0){
      /*fprintf(stderr,"Invalid bin #%d\n\r",b); */
      continue;
    }

    bin[b].theta /= bin[b].dtheta;
    phi0   /= bin[b].offset;
    phi0   += ((float)shift + (float)p0) * 360./32.;
    if(phi0 < 0.)
      phi0 +=360.;
    if(phi0 >= 360.)
      phi0 -=360.;
    bin[b].dtheta /= num_phi[b];
    bin[b].dphi   /= num_theta[b];
    bin[b].geom   /= num_phi[b];    
    /* build up phi across energy bins */
    for (e=0; e<bin[b].ne; e++)
      bin[b].phi[e] = phi0 +	(.5 - (e+.5)/bin[b].ne)*dphi;
  }

  offset = 0;
  for(b=0;b<nbins;b++){
    bin[b].offset = offset;
    offset += bin[b].ne;
  }
  if(offset != nsamples)
    fprintf(stderr,"3d Initialization error, offset: %d, nsamples: %d\n\r",
	    offset, nsamples);
  
  memcpy(&prev_mapinfo,&mapinfo,sizeof(mapinfo));
  return(1);
}

int decom_map_ph97(packet *pk, idl_ph97 *map){
  int b,e,i,p,t;
  int seq,offset,ndata,nsamples,npackets,max_bytes_per_packet=390;
  static int status, double_sweep, num_phi[MAX_PH_BINS];
  static float alldata[1530];
  double dt,nrg[15],dnrg[15];
  uint    ps=0;
  static  uint ptmap[32*32];
  static BinRange bin[97]; 
  uchar *data,shift,p0=28;
  float flux,n0,phi;
  packet dpk;

  i=offset=seq=b=e=p=t=flux=ndata=npackets=0;
  nsamples=1530;

  decompress_burst_packet(&dpk,pk);
  pk = &dpk;

  if(pk->spin != map->spin || status<=0){
    /* This is executed only on the first of "npackets" calls.
       It gets all of the angular and energy information, and 
       sets things such as the time and spin. */
    PCFG *pcfg;   

    for (i=0;i<nsamples;i++) alldata[i]=0.;

    pcfg = get_PCFG(pk->time);

    if((pcfg->norm_cfg.esa_swp_select==0x7780) &&    
       (pcfg->norm_cfg.select_sector ==0x0343)  &&
       (pcfg->norm_cfg.esa_pha_basech==0x0040)){
      double_sweep=1;
      /* this marks the weird double sweep of the top half energy channels */
    } else {
      double_sweep=0;
    }

#ifdef PH_DEBUG
    printf("esa_swp_select: %04x\n\r",pcfg->norm_cfg.esa_swp_select);
    printf("select_sector : %04x\n\r",pcfg->norm_cfg.select_sector );
    printf("esa_pha_basech: %04x\n\r",pcfg->norm_cfg.esa_pha_basech);
    printf("bts_c_val     : %04x\n\r",pcfg->norm_cfg.bts_c_val);
    printf("double_sweep  : %4d\n\r",double_sweep);
#endif

    map->time        = pk->time;
    /* get_3d_integ_t(pk): number of summed spins */
    map->integ_t     = get_3d_integ_t(pk) * pcfg->spin_period;
    map->delta_t     = get_3d_delta_t(pk) * pcfg->spin_period;
    map->end_time    = map->time + map->integ_t;
    map->geom_factor = pcfg->ph_geom;
    map->double_sweep = (uchar) double_sweep;
    map->spin        = pk->spin;
    /* these are set in idl
    map->mass        = 1836*5.6856591e-6;
    map->nbins       = 97;
    map->nenergy     = 15; */

    dt       = map->integ_t / 32 / 16; /* nspins*spin period/nphis/nenergies */
    shift    = (pk->instseq >> 3) & 0x1f;
    map->shift = shift;
    npackets = 4;   
    status   = npackets;  /*decrement variable, counts packets left to gather*/
      
    for(i=0;i<30*4;i++){
      map->dac_code[i] = pcfg->phdac_tbl[i];
      map->volts[i]    = pcfg->ph_volts_tbl[i]; 
    }
    
    make_ph_ptmap(pcfg,shift,map->mapcode,num_phi,ptmap,bin);

    for(b=0;b<map->nbins;b++){  /* initalize data structure */
      map->domega[b] = 0.;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = 0.;
	map->dnrg [b][e] = 0.;
	map->theta[b][e] = 0.;
	map->phi  [b][e] = 0.;
	map->dtheta[b][e] = 0.;
        map->dphi  [b][e] = 0.;
	map->flux [b][e] = 0.;
	map->gf   [b][e] = 0.;
	map->dt   [b][e] = 0.;
      }
    }
    
    get_esteps_ph(nrg,map->nenergy,MIDDLE,pcfg);
    get_esteps_ph(dnrg,map->nenergy,WIDTH,pcfg);

  }

  /* this is the *start* of the section executed on each call */
  seq    = pk->instseq & 0x07;
  ndata  = max_bytes_per_packet;
  offset = seq * ndata;
  if (ndata>(int)pk->dsize) ndata = pk->dsize;
  data   = pk->data;

  for(i=0;i<ndata;i++){
    /*if(offset >= nsamples){ */
      if(offset >= (nsamples + 5))  {   /*The previous line should work */
      char instStr[12];
      strcpy (instStr, "PESA HIGH");
      fprintf(stderr,"3D binning error: offset(%d) >= nsamples(%d)\n\r",
	      offset, nsamples);
      fprintf(stderr,"    ndata: %d, map type: %s\n\r", ndata, instStr);
    } else {
      if(offset < nsamples)
	alldata[offset++] = decomp19_8(*data++);
    }
  }
  status--;
  /* this is the **end** of the section executed on each call */
  
  if (status == 0) {   /* if all packets received, fill the data structure */

    uint ts=0;
    for(t=0;t<32;t++){ /* fill the phi-theta map */
      ts = tsect[t];
      for(p=0;p<32;p++){
        if(t<24)
        	map->pt_map[t*32+p] = ptmap[t*32+p];
	b = ptmap[ts*32+p];
	ps = (p - p0) & 0x1f;
	if (b > 97) 
	  printf("invalid map entry: t: %d, p: %d, ts: %d, ps: %d, t*32+ps: %d, b: %d\n\r",
		 t,p,ts,ps,t*32+ps,b);
      }
    }
    
    for(b=0;b<map->nbins;b++){
      offset         = bin[b].offset;
      map->domega[b] = bin[b].domega;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = nrg[e];
	map->dnrg [b][e] = dnrg[e];
	map->theta[b][e] = bin[b].theta;
	map->phi  [b][e] = bin[b].phi[e];
        map->dtheta[b][e] = bin[b].dtheta;
        map->dphi  [b][e] = bin[b].dphi;
	map->gf   [b][e] = bin[b].geom;
	map->dt   [b][e] = num_phi[b] * dt;
	map->dvolume[b][e] = dnrg[e]/nrg[e]*map->domega[b];
	if (map->nenergy==bin[b].ne) {
	  flux  = alldata[offset+e];
	  phi 	= bin[b].phi[e];
	} else if (map->nenergy==bin[b].ne/2) {
	  flux  = alldata[offset+2*e];
	  flux += alldata[offset+2*e+1];
	  phi	= bin[b].phi[2*e];
	  phi  += bin[b].phi[2*e+1];
	  phi	= phi/2.;
	} else {
	  flux  = 0;
	  /*printf("bin: %d: energy steps=%d, should be %d or %d\n\r",
		 b,bin[b].ne,map->nenergy,2*map->nenergy);*/
	}
	map->flux[b][e]  = flux;
	map->phi[b][e] = phi;
      }
    }

    if(double_sweep){
      int dsbins[4];
      int midbin;                       /* ramping bin between double sweeps */

      midbin=(map->nenergy-2+map->nenergy%2)/2;

      for(t=0;t<16;t++){
        dsbins[0] = ptmap[t*32+26];
        dsbins[1] = ptmap[t*32+27];
        dsbins[2] = ptmap[t*32+28];
        dsbins[3] = ptmap[t*32+29];
        if((dsbins[0] == dsbins[1]) || (dsbins[1] == dsbins[3]))
          for(i=0;i<6;i++)
	  for(e=midbin;e<map->nenergy;e++) {
	    map->flux[dsbins[1]][e] = NaN;        /* blank corrupted bins */
	    map->bins[dsbins[1]][e] = 0;
	  }
	else
	  for(i=0;i<16;i++){
	    map->nrg[dsbins[1]][midbin] = NaN;
	    map->dnrg[dsbins[1]][midbin] = NaN;
	    map->bins[dsbins[1]][midbin] = 0;
	    for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	      map->nrg [dsbins[1]][e] = map->nrg [dsbins[1]][e-midbin-1]; /*fix  e */
	      map->dnrg[dsbins[1]][e] = map->dnrg[dsbins[1]][e-midbin-1]; /*fix de */
	    }
	  }
        if(dsbins[1] != dsbins[2]) {
          if((dsbins[0] == dsbins[2]) || (dsbins[2] == dsbins[3]))
            for(i=0;i<6;i++)
  	    for(e=midbin;e<map->nenergy;e++) {
	      map->flux[dsbins[2]][e] = NaN;        /* blank corrupted bins */
	      map->bins[dsbins[2]][e] = 0;
	    }
	  else
	    for(i=0;i<16;i++){
	      map->nrg[dsbins[2]][midbin] = NaN;
	      map->dnrg[dsbins[2]][midbin] = NaN;
	      map->bins[dsbins[2]][midbin] = 0;
	      for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	        map->nrg [dsbins[2]][e] = map->nrg [dsbins[2]][e-midbin-1]; /*fix  e */
	        map->dnrg[dsbins[2]][e] = map->dnrg[dsbins[2]][e-midbin-1]; /*fix de */
	      }
	    }
	}
      }
    }
  
    map->valid=1;
  }

  return(1);  
}

int decom_map_ph121(packet *pk, idl_ph121 *map){
  int     b,e,i,p,t;
  int     seq,offset,ndata,nsamples,npackets,max_bytes_per_packet=390;
  static  int      status, double_sweep, num_phi[MAX_PH_BINS];
  static  float    alldata[1890];
  double  dt,nrg[15],dnrg[15];
  uint    ps=0;
  static  uint ptmap[32*32];
  static  BinRange bin[121]; 
  uchar  *data,shift,p0=28;
  float   flux,n0,phi;
  packet  dpk;

  i=offset=seq=b=e=p=t=flux=ndata=npackets=0;
  nsamples=1890;

  decompress_burst_packet(&dpk,pk);
  pk = &dpk;

  if(pk->spin != map->spin || status<=0){
    /* This is executed only on the first of "npackets" calls.
       It gets all of the angular and energy information, and 
       sets things such as the time and spin. */
    PCFG *pcfg;   

    for (i=0;i<nsamples;i++) alldata[i]=0.;

    pcfg = get_PCFG(pk->time);
  
    if((pcfg->norm_cfg.esa_swp_select==0x7780) &&    
       (pcfg->norm_cfg.select_sector ==0x0343)  &&
       (pcfg->norm_cfg.esa_pha_basech==0x0040)){
      double_sweep=1;
      /* this marks the weird double sweep of the top half energy channels */
    } else {
      double_sweep=0;
    }

#ifdef PH_DEBUG
    printf("esa_swp_select: %04x\n\r",pcfg->norm_cfg.esa_swp_select);
    printf("select_sector : %04x\n\r",pcfg->norm_cfg.select_sector );
    printf("esa_pha_basech: %04x\n\r",pcfg->norm_cfg.esa_pha_basech);
    printf("bts_c_val     : %04x\n\r",pcfg->norm_cfg.bts_c_val);
    printf("double_sweep  : %4d\n\r",double_sweep);
#endif 
  
    map->time        = pk->time;
    /* get_3d_integ_t(pk): number of summed spins */
    map->integ_t     = get_3d_integ_t(pk) * pcfg->spin_period;
    map->delta_t     = get_3d_delta_t(pk) * pcfg->spin_period;
    map->end_time    = map->time + map->integ_t;
    map->geom_factor = pcfg->ph_geom;
    map->double_sweep = (uchar) double_sweep;
    map->spin        = pk->spin;
    /* these are set in idl
    map->mass        = 1836*5.6856591e-6;
    map->nbins       = 121;
    map->nenergy     = 15; */

    dt       = map->integ_t / 32 / 16; /* nspins*spin period/nphis/nenergies */
    shift    = (pk->instseq >> 3) & 0x1f;
    map->shift = shift;
    npackets = 5;   
    status   = npackets; /*decrement variable, counts packets left to gather */
    
    for(i=0;i<30*4;i++){
      map->dac_code[i] = pcfg->phdac_tbl[i];
      map->volts[i]    = pcfg->ph_volts_tbl[i]; 
    }
    
    make_ph_ptmap(pcfg,shift,map->mapcode,num_phi,ptmap,bin);

    for(b=0;b<map->nbins;b++){  /* initalize data structure */ 
      map->domega[b] = 0.;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = 0.;
	map->dnrg [b][e] = 0.;
	map->theta[b][e] = 0.;
	map->phi  [b][e] = 0.;
        map->dtheta[b][e] = 0.;
        map->dphi  [b][e] = 0.;
	map->flux [b][e] = 0.;  /* Nan ? */
	map->gf   [b][e] = 0.;
	map->dt   [b][e] = 0.;
      }
    }

    get_esteps_ph(nrg,map->nenergy,MIDDLE,pcfg);
    get_esteps_ph(dnrg,map->nenergy,WIDTH,pcfg);

  }

  /* this is the *start* of the section executed on each call */
  seq    = pk->instseq & 0x07;
  ndata  = max_bytes_per_packet;
  offset = seq * ndata;
  if (ndata>(int)pk->dsize) ndata = pk->dsize;
  data   = pk->data;

  for(i=0;i<ndata;i++){
    /*if(offset >= nsamples){ */
      if(offset >= (nsamples + 5))  {   /* The previous line should work */
      char instStr[12];
      strcpy (instStr, "PESA HIGH");
      fprintf(stderr,"3D binning error: offset(%d) >= nsamples(%d)\n\r",
	      offset, nsamples);
      fprintf(stderr,"    ndata: %d, map type: %s\n\r", ndata, instStr);
    } else {
      if(offset < nsamples)
	alldata[offset++] = decomp19_8(*data++);
    }
  }
  status--;
  /* this is the **end** of the section executed on each call */
  
  if (status == 0) {   /* if all packets received, fill the data structure */

    uint ts=0;
    for(t=0;t<32;t++){ /* fill the phi-theta map */
      ts = tsect[t];
      for(p=0;p<32;p++){
        if(t<24)
        	map->pt_map[t*32+p] = ptmap[t*32+p];
	b = ptmap[ts*32+p];
	ps = (p - p0) & 0x1f;
      }
    }
    
    for(b=0;b<map->nbins;b++){
      offset         = bin[b].offset;
      map->domega[b] = bin[b].domega;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = nrg[e];
	map->dnrg [b][e] = dnrg[e];
	map->theta[b][e] = bin[b].theta;
        map->dphi  [b][e] = bin[b].dphi;
        map->dtheta[b][e] = bin[b].dtheta;
	map->gf   [b][e] = bin[b].geom;
	if(b==120)
		map->dt[b][e] = 8 * dt;
	else
		map->dt[b][e] = num_phi[b] * dt;
	map->dvolume[b][e] = dnrg[e]/nrg[e]*map->domega[b];
	if (map->nenergy==bin[b].ne) {
	  flux  = alldata[offset+e];  
	  phi 	= bin[b].phi[e];
	} else if (map->nenergy==bin[b].ne/2) {
	  flux  = alldata[offset+2*e];
	  flux += alldata[offset+2*e+1];
	  phi	= bin[b].phi[2*e];
	  phi  += bin[b].phi[2*e+1];
	  phi	= phi/2.;
	} else {
	  flux  = 0;
	  /*printf("bin: %d: energy steps=%d, should be %d or %d\n\r",
		 b,bin[b].ne,map->nenergy,2*map->nenergy);*/
	}
	map->flux[b][e]  = flux;
 	map->phi[b][e] = phi;
     }
    }

    if(double_sweep){
      int dsbins[4];
      int midbin;                       /* ramping bin between double sweeps */

      midbin=(map->nenergy-2+map->nenergy%2)/2;

      for(t=0;t<16;t++){
        dsbins[0] = ptmap[t*32+26];
        dsbins[1] = ptmap[t*32+27];
        dsbins[2] = ptmap[t*32+28];
        dsbins[3] = ptmap[t*32+29];
        if((dsbins[0] == dsbins[1]) || (dsbins[1] == dsbins[3]))
          for(i=0;i<6;i++)
	  for(e=midbin;e<map->nenergy;e++) {
	    map->flux[dsbins[1]][e] = NaN;        /* blank corrupted bins */
	    map->bins[dsbins[1]][e] = 0;        /* blank corrupted bins */
	  }
	else
	  for(i=0;i<16;i++){
	    map->nrg[dsbins[1]][midbin] = NaN;
	    map->dnrg[dsbins[1]][midbin] = NaN;
	    map->bins[dsbins[1]][midbin] = 0;        /* blank corrupted bins */
	    for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	      map->nrg [dsbins[1]][e] = map->nrg [dsbins[1]][e-midbin-1]; /*fix  e */
	      map->dnrg[dsbins[1]][e] = map->dnrg[dsbins[1]][e-midbin-1]; /*fix de */
	    }
	  }
        if(dsbins[1] != dsbins[2]) {
          if((dsbins[0] == dsbins[2]) || (dsbins[2] == dsbins[3]))
            for(i=0;i<6;i++)
  	    for(e=midbin;e<map->nenergy;e++) {
	      map->flux[dsbins[2]][e] = NaN;        /* blank corrupted bins */
	      map->bins[dsbins[2]][e] = 0;        /* blank corrupted bins */
	    }
	  else
	    for(i=0;i<16;i++){
	      map->nrg[dsbins[2]][midbin] = NaN;
	      map->dnrg[dsbins[2]][midbin] = NaN;
	      map->bins[dsbins[2]][midbin] = 0;        /* blank corrupted bins */
	      for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	        map->nrg [dsbins[2]][e] = map->nrg [dsbins[2]][e-midbin-1]; /*fix  e */
	        map->dnrg[dsbins[2]][e] = map->dnrg[dsbins[2]][e-midbin-1]; /*fix de */
	      }
	    }
	}
      }
    }

    map->valid=1;
  }

  return(1);  
}

int decom_map_ph56(packet *pk, idl_ph56 *map){
  int b,e,i,p,t;
  int seq,offset,ndata,nsamples,npackets,max_bytes_per_packet=390;
  static int status, double_sweep, num_phi[MAX_PH_BINS];
  static float alldata[960];
  double dt,nrg[15],dnrg[15];
  uint    ps=0;
  static  uint ptmap[32*32];
  static BinRange bin[56]; 
  uchar *data,shift,p0=28;
  float flux,n0,phi;
  packet dpk;

  i=offset=seq=b=e=p=t=flux=ndata=npackets=0;
  nsamples=960;

  decompress_burst_packet(&dpk,pk);
  pk = &dpk;

  if(pk->spin != map->spin || status<=0){
    /* This is executed only on the first of "npackets" calls.
       It gets all of the angular and energy information, and 
       sets things such as the time and spin. */
    PCFG *pcfg;   

    for (i=0;i<nsamples;i++) alldata[i]=0.;

    pcfg = get_PCFG(pk->time);

    if((pcfg->norm_cfg.esa_swp_select==0x7780) &&    
       (pcfg->norm_cfg.select_sector ==0x0343)  &&
       (pcfg->norm_cfg.esa_pha_basech==0x0040)){
      double_sweep=1;
      /* this marks the weird double sweep of the top half energy channels */
    } else {
      double_sweep=0;
    }

#ifdef PH_DEBUG
    printf("esa_swp_select: %04x\n\r",pcfg->norm_cfg.esa_swp_select);
    printf("select_sector : %04x\n\r",pcfg->norm_cfg.select_sector );
    printf("esa_pha_basech: %04x\n\r",pcfg->norm_cfg.esa_pha_basech);
    printf("bts_c_val     : %04x\n\r",pcfg->norm_cfg.bts_c_val);
    printf("double_sweep  : %4d\n\r",double_sweep);
#endif

    map->time        = pk->time;
    /* get_3d_integ_t(pk): number of summed spins */
    map->integ_t     = get_3d_integ_t(pk) * pcfg->spin_period;
    map->delta_t     = get_3d_delta_t(pk) * pcfg->spin_period;
    map->end_time    = map->time + map->integ_t;
    map->geom_factor = pcfg->ph_geom;
    map->double_sweep = (uchar) double_sweep;
    map->spin        = pk->spin;
    /* these are set in idl
    map->mass        = 1836*5.6856591e-6;
    map->nbins       = 56;
    map->nenergy     = 15; */

    dt       = map->integ_t / 32 / 16; /* nspins*spin period/nphis/nenergies */
    shift    = (pk->instseq >> 3) & 0x1f;
    map->shift = shift;
    npackets = 3;   
    status   = npackets;  /*decrement variable, counts packets left to gather*/
      
    for(i=0;i<30*4;i++){
      map->dac_code[i] = pcfg->phdac_tbl[i];
      map->volts[i]    = pcfg->ph_volts_tbl[i]; 
    }
    
    make_ph_ptmap(pcfg,shift,map->mapcode,num_phi,ptmap,bin);

    for(b=0;b<map->nbins;b++){  /* initalize data structure */
      map->domega[b] = 0.;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = 0.;
	map->dnrg [b][e] = 0.;
	map->theta[b][e] = 0.;
	map->phi  [b][e] = 0.;
        map->dtheta[b][e] = 0.;
        map->dphi  [b][e] = 0.;
	map->flux [b][e] = 0.;
	map->gf   [b][e] = 0.;
	map->dt   [b][e] = 0.;
      }
    }
    
    get_esteps_ph(nrg,map->nenergy,MIDDLE,pcfg);
    get_esteps_ph(dnrg,map->nenergy,WIDTH,pcfg);

  }

  /* this is the *start* of the section executed on each call */
  seq    = pk->instseq & 0x07;
  ndata  = max_bytes_per_packet;
  offset = seq * ndata;
  if (ndata>(int)pk->dsize) ndata = pk->dsize;
  data   = pk->data;

  for(i=0;i<ndata;i++){
    /*if(offset >= nsamples){*/
      if(offset >= (nsamples + 5))  {   /* The previous line should work */
      char instStr[12];
      strcpy (instStr, "PESA HIGH");
      fprintf(stderr,"3D binning error: offset(%d) >= nsamples(%d)\n\r",
	      offset, nsamples);
      fprintf(stderr,"    ndata: %d, map type: %s\n\r", ndata, instStr);
    } else {
      if(offset < nsamples)
	alldata[offset++] = decomp19_8(*data++);
    }
  }
  status--;
  /* this is the **end** of the section executed on each call */
  
  if (status == 0) {   /* if all packets received, fill the data structure */

    uint ts=0;
    for(t=0;t<32;t++){ /* fill the phi-theta map */
      ts = tsect[t];
      for(p=0;p<32;p++){
        if(t<24)
        	map->pt_map[t*32+p] = ptmap[t*32+p];
	b = ptmap[ts*32+p];
	ps = (p - p0) & 0x1f;
	if (b > 56) 
	  printf("invalid map entry: t: %d, p: %d, ts: %d, ps: %d, t*32+ps: %d, b: %d\n\r",
		 t,p,ts,ps,t*32+ps,b);
      }
    }
    
    for(b=0;b<map->nbins;b++){
      offset         = bin[b].offset;
      map->domega[b] = bin[b].domega;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = nrg[e];
	map->dnrg [b][e] = dnrg[e];
	map->theta[b][e] = bin[b].theta;
        map->dtheta[b][e] = bin[b].dtheta;
        map->dphi  [b][e] = bin[b].dphi;
	map->gf   [b][e] = bin[b].geom;
	map->dt   [b][e] = num_phi[b] * dt;
	map->dvolume[b][e] = dnrg[e]/nrg[e]*map->domega[b];
	if (map->nenergy==bin[b].ne) {
	  flux  = alldata[offset+e];  
	  phi 	= bin[b].phi[e];
	} else if (map->nenergy==bin[b].ne/2) {
	  flux  = alldata[offset+2*e];
	  flux += alldata[offset+2*e+1];
	  phi	= bin[b].phi[2*e];
	  phi  += bin[b].phi[2*e+1];
	  phi	= phi/2.;
	} else {
	  flux  = 0;
	  /*printf("bin: %d: energy steps=%d, should be %d or %d\n\r",
		 b,bin[b].ne,map->nenergy,2*map->nenergy);*/
	}
	map->flux[b][e]  = flux;
	map->phi[b][e] = phi;
      }
    }

    if(double_sweep){
      int dsbins[4];
      int midbin;                       /* ramping bin between double sweeps */

      midbin=(map->nenergy-2+map->nenergy%2)/2;

      for(t=0;t<16;t++){
        dsbins[0] = ptmap[t*32+26];
        dsbins[1] = ptmap[t*32+27];
        dsbins[2] = ptmap[t*32+28];
        dsbins[3] = ptmap[t*32+29];
        if((dsbins[0] == dsbins[1]) || (dsbins[1] == dsbins[3]))
          for(i=0;i<6;i++)
	  for(e=midbin;e<map->nenergy;e++) {
	    map->flux[dsbins[1]][e] = NaN;        /* blank corrupted bins */
	    map->bins[dsbins[1]][e] = 0;        /* blank corrupted bins */
	  }
	else
	  for(i=0;i<16;i++){
	    map->nrg[dsbins[1]][midbin] = NaN;
	    map->dnrg[dsbins[1]][midbin] = NaN;
	    map->bins[dsbins[1]][midbin] = 0;        /* blank corrupted bins */
	    for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	      map->nrg [dsbins[1]][e] = map->nrg [dsbins[1]][e-midbin-1]; /*fix  e */
	      map->dnrg[dsbins[1]][e] = map->dnrg[dsbins[1]][e-midbin-1]; /*fix de */
	    }
	  }
        if(dsbins[1] != dsbins[2]) {
          if((dsbins[0] == dsbins[2]) || (dsbins[2] == dsbins[3]))
            for(i=0;i<6;i++)
  	    for(e=midbin;e<map->nenergy;e++) {
	      map->flux[dsbins[2]][e] = NaN;        /* blank corrupted bins */
	      map->bins[dsbins[2]][e] = 0;        /* blank corrupted bins */
	    }
	  else
	    for(i=0;i<16;i++){
	      map->nrg[dsbins[2]][midbin] = NaN;
	      map->dnrg[dsbins[2]][midbin] = NaN;
	      map->bins[dsbins[2]][midbin] = 0;        /* blank corrupted bins */
	      for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	        map->nrg [dsbins[2]][e] = map->nrg [dsbins[2]][e-midbin-1]; /*fix  e */
	        map->dnrg[dsbins[2]][e] = map->dnrg[dsbins[2]][e-midbin-1]; /*fix de */
	      }
	    }
	}
      }
    }

    map->valid=1;
  }

  return(1);  
}

int decom_map_ph65(packet *pk, idl_ph65 *map){
  int b,e,i,p,t;
  int seq,offset,ndata,nsamples,npackets,max_bytes_per_packet=390;
  static int status, double_sweep, num_phi[MAX_PH_BINS];
  static float alldata[1050];
  double dt,nrg[15],dnrg[15];
  uint    ps=0;
  static  uint ptmap[32*32];
  static BinRange bin[65]; 
  uchar *data,shift,p0=28;
  float flux,n0,phi;
  packet dpk;

  i=offset=seq=b=e=p=t=flux=ndata=npackets=0;
  nsamples=1050;

  decompress_burst_packet(&dpk,pk);
  pk = &dpk;

  if(pk->spin != map->spin || status<=0){
    /* This is executed only on the first of "npackets" calls.
       It gets all of the angular and energy information, and 
       sets things such as the time and spin. */
    PCFG *pcfg;   

    for (i=0;i<nsamples;i++) alldata[i]=0.;

    pcfg = get_PCFG(pk->time);

    if((pcfg->norm_cfg.esa_swp_select==0x7780) &&    
       (pcfg->norm_cfg.select_sector ==0x0343)  &&
       (pcfg->norm_cfg.esa_pha_basech==0x0040)){
      double_sweep=1;
      /* this marks the weird double sweep of the top half energy channels */
    } else {
      double_sweep=0;
    }

#ifdef PH_DEBUG
    printf("esa_swp_select: %04x\n\r",pcfg->norm_cfg.esa_swp_select);
    printf("select_sector : %04x\n\r",pcfg->norm_cfg.select_sector );
    printf("esa_pha_basech: %04x\n\r",pcfg->norm_cfg.esa_pha_basech);
    printf("bts_c_val     : %04x\n\r",pcfg->norm_cfg.bts_c_val);
    printf("double_sweep  : %4d\n\r",double_sweep);
#endif

    map->time        = pk->time;
    /* get_3d_integ_t(pk): number of summed spins */
    map->integ_t     = get_3d_integ_t(pk) * pcfg->spin_period;
    map->delta_t     = get_3d_delta_t(pk) * pcfg->spin_period;
    map->end_time    = map->time + map->integ_t;
    map->geom_factor = pcfg->ph_geom;
    map->double_sweep = (uchar) double_sweep;
    map->spin        = pk->spin;
    /* these are set in idl
    map->mass        = 1836*5.6856591e-6;
    map->nbins       = 65;
    map->nenergy     = 15; */

    dt       = map->integ_t / 32 / 16; /* nspins*spin period/nphis/nenergies */
    shift    = (pk->instseq >> 3) & 0x1f;
    map->shift = shift;
    npackets = 3;   
    status   = npackets;  /*decrement variable, counts packets left to gather*/
      
    for(i=0;i<30*4;i++){
      map->dac_code[i] = pcfg->phdac_tbl[i];
      map->volts[i]    = pcfg->ph_volts_tbl[i]; 
    }
    
    make_ph_ptmap(pcfg,shift,map->mapcode,num_phi,ptmap,bin);

    for(b=0;b<map->nbins;b++){  /* initalize data structure */
      map->domega[b] = 0.;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = 0.;
	map->dnrg [b][e] = 0.;
	map->theta[b][e] = 0.;
	map->phi  [b][e] = 0.;
        map->dtheta[b][e] = 0.;
        map->dphi  [b][e] = 0.;
	map->flux [b][e] = 0.;
	map->gf   [b][e] = 0.;
	map->dt   [b][e] = 0.;
      }
    }
    
    get_esteps_ph(nrg,map->nenergy,MIDDLE,pcfg);
    get_esteps_ph(dnrg,map->nenergy,WIDTH,pcfg);

  }

  /* this is the *start* of the section executed on each call */
  seq    = pk->instseq & 0x07;
  ndata  = max_bytes_per_packet;
  offset = seq * ndata;
  if (ndata>(int)pk->dsize) ndata = pk->dsize;
  data   = pk->data;

  for(i=0;i<ndata;i++){
    /*if(offset >= nsamples){ */
      if(offset >= (nsamples + 5))  {   /* The previous line should work */
      char instStr[12];
      strcpy (instStr, "PESA HIGH");
      fprintf(stderr,"3D binning error: offset(%d) >= nsamples(%d)\n\r",
	      offset, nsamples);
      fprintf(stderr,"    ndata: %d, map type: %s\n\r", ndata, instStr);
    } else {
      if(offset < nsamples)
	alldata[offset++] = decomp19_8(*data++);
    }
  }
  status--;
  /* this is the **end** of the section executed on each call */
  
  if (status == 0) {   /* if all packets received, fill the data structure */

    uint ts=0;
    for(t=0;t<32;t++){ /* fill the phi-theta map */
      ts = tsect[t];
      for(p=0;p<32;p++){
        if(t<24)
        	map->pt_map[t*32+p] = ptmap[t*32+p];
	b = ptmap[ts*32+p];
	ps = (p - p0) & 0x1f;
	if (b > 65) 
	  printf("invalid map entry: t: %d, p: %d, ts: %d, ps: %d, t*32+ps: %d, b: %d\n\r",
		 t,p,ts,ps,t*32+ps,b);
      }
    }
    
    for(b=0;b<map->nbins;b++){
      offset         = bin[b].offset;
      map->domega[b] = bin[b].domega;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = nrg[e];
	map->dnrg [b][e] = dnrg[e];
	map->theta[b][e] = bin[b].theta;
        map->dtheta[b][e] = bin[b].dtheta;
        map->dphi  [b][e] = bin[b].dphi;
	map->gf   [b][e] = bin[b].geom;
	map->dt   [b][e] = num_phi[b] * dt;
	map->dvolume[b][e] = dnrg[e]/nrg[e]*map->domega[b];
	if (map->nenergy==bin[b].ne) {
	  flux  = alldata[offset+e];  
	  phi 	= bin[b].phi[e];
	} else if (map->nenergy==bin[b].ne/2) {
	  flux  = alldata[offset+2*e];
	  flux += alldata[offset+2*e+1];
	  phi	= bin[b].phi[2*e];
	  phi  += bin[b].phi[2*e+1];
	  phi	= phi/2.;
	} else {
	  flux  = 0;
	  /*printf("bin: %d: energy steps=%d, should be %d or %d\n\r",
		 b,bin[b].ne,map->nenergy,2*map->nenergy);*/
	}
	map->flux[b][e]  = flux;
	map->phi[b][e] = phi;
      }
    }

    if(double_sweep){
      int dsbins[4];
      int midbin;                       /* ramping bin between double sweeps */

      midbin=(map->nenergy-2+map->nenergy%2)/2;

      for(t=0;t<16;t++){
        dsbins[0] = ptmap[t*32+26];
        dsbins[1] = ptmap[t*32+27];
        dsbins[2] = ptmap[t*32+28];
        dsbins[3] = ptmap[t*32+29];
        if((dsbins[0] == dsbins[1]) || (dsbins[1] == dsbins[3]))
          for(i=0;i<6;i++)
	  for(e=midbin;e<map->nenergy;e++) {
	    map->flux[dsbins[1]][e] = NaN;        /* blank corrupted bins */
	    map->bins[dsbins[1]][e] = 0;        /* blank corrupted bins */
	  }
	else
	  for(i=0;i<16;i++){
	    map->nrg[dsbins[1]][midbin] = NaN;
	    map->dnrg[dsbins[1]][midbin] = NaN;
	    map->bins[dsbins[1]][midbin] = 0;        /* blank corrupted bins */
	    for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	      map->nrg [dsbins[1]][e] = map->nrg [dsbins[1]][e-midbin-1]; /*fix  e */
	      map->dnrg[dsbins[1]][e] = map->dnrg[dsbins[1]][e-midbin-1]; /*fix de */
	    }
	  }
        if(dsbins[1] != dsbins[2]) {
          if((dsbins[0] == dsbins[2]) || (dsbins[2] == dsbins[3]))
            for(i=0;i<6;i++)
  	    for(e=midbin;e<map->nenergy;e++) {
	      map->flux[dsbins[2]][e] = NaN;        /* blank corrupted bins */
	      map->bins[dsbins[2]][e] = 0;        /* blank corrupted bins */
	    }
	  else
	    for(i=0;i<16;i++){
	      map->nrg[dsbins[2]][midbin] = NaN;
	      map->dnrg[dsbins[2]][midbin] = NaN;
	      map->bins[dsbins[2]][midbin] = 0;        /* blank corrupted bins */
	      for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	        map->nrg [dsbins[2]][e] = map->nrg [dsbins[2]][e-midbin-1]; /*fix  e */
	        map->dnrg[dsbins[2]][e] = map->dnrg[dsbins[2]][e-midbin-1]; /*fix de */
	      }
	    }
	}
      }
    }

    map->valid=1;
  }

  return(1);  
}



int decom_map_ph88(packet *pk, idl_ph88 *map){
  int b,e,i,p,t;
  int seq,offset,ndata,nsamples,npackets,max_bytes_per_packet=390;
  static int status, double_sweep, num_phi[MAX_PH_BINS];
  static float alldata[1380];
  double dt,nrg[15],dnrg[15];
  uint    ps=0;
  static  uint ptmap[32*32];
  static BinRange bin[88]; 
  uchar *data,shift,p0=28;
  float flux,n0,phi;
  packet dpk;

  i=offset=seq=b=e=p=t=flux=ndata=npackets=0;
  nsamples=1380;		/* this and alldata may be 1470 */

  decompress_burst_packet(&dpk,pk);
  pk = &dpk;

  if(pk->spin != map->spin || status<=0){
    /* This is executed only on the first of "npackets" calls.
       It gets all of the angular and energy information, and 
       sets things such as the time and spin. */
    PCFG *pcfg;   

    for (i=0;i<nsamples;i++) alldata[i]=0.;

    pcfg = get_PCFG(pk->time);

    if((pcfg->norm_cfg.esa_swp_select==0x7780) &&    
       (pcfg->norm_cfg.select_sector ==0x0343)  &&
       (pcfg->norm_cfg.esa_pha_basech==0x0040)){
      double_sweep=1;
      /* this marks the weird double sweep of the top half energy channels */
    } else {
      double_sweep=0;
    }

#ifdef PH_DEBUG
    printf("esa_swp_select: %04x\n\r",pcfg->norm_cfg.esa_swp_select);
    printf("select_sector : %04x\n\r",pcfg->norm_cfg.select_sector );
    printf("esa_pha_basech: %04x\n\r",pcfg->norm_cfg.esa_pha_basech);
    printf("bts_c_val     : %04x\n\r",pcfg->norm_cfg.bts_c_val);
    printf("double_sweep  : %4d\n\r",double_sweep);
#endif

    map->time        = pk->time;
    /* get_3d_integ_t(pk): number of summed spins */
    map->integ_t     = get_3d_integ_t(pk) * pcfg->spin_period;
    map->delta_t     = get_3d_delta_t(pk) * pcfg->spin_period;
    map->end_time    = map->time + map->integ_t;
    map->geom_factor = pcfg->ph_geom;
    map->double_sweep = (uchar) double_sweep;
    map->spin        = pk->spin;
    /* these are set in idl
    map->mass        = 1836*5.6856591e-6;
    map->nbins       = 88;
    map->nenergy     = 15; */

    dt       = map->integ_t / 32 / 16; /* nspins*spin period/nphis/nenergies */
    shift    = (pk->instseq >> 3) & 0x1f;
    map->shift = shift;
    npackets = 4;   
    status   = npackets;  /*decrement variable, counts packets left to gather*/
      
    for(i=0;i<30*4;i++){
      map->dac_code[i] = pcfg->phdac_tbl[i];
      map->volts[i]    = pcfg->ph_volts_tbl[i]; 
    }
    
    make_ph_ptmap(pcfg,shift,map->mapcode,num_phi,ptmap,bin);

    for(b=0;b<map->nbins;b++){  /* initalize data structure */
      map->domega[b] = 0.;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = 0.;
	map->dnrg [b][e] = 0.;
	map->theta[b][e] = 0.;
	map->phi  [b][e] = 0.;
        map->dtheta[b][e] = 0.;
        map->dphi  [b][e] = 0.;
	map->flux [b][e] = 0.;
	map->gf   [b][e] = 0.;
	map->dt   [b][e] = 0.;
      }
    }
    
    get_esteps_ph(nrg,map->nenergy,MIDDLE,pcfg);
    get_esteps_ph(dnrg,map->nenergy,WIDTH,pcfg);

  }

  /* this is the *start* of the section executed on each call */
  seq    = pk->instseq & 0x07;
  ndata  = max_bytes_per_packet;
  offset = seq * ndata;
  if (ndata>(int)pk->dsize) ndata = pk->dsize;
  data   = pk->data;

  for(i=0;i<ndata;i++){
    /*if(offset >= nsamples){ */
      if(offset >= (nsamples + 5))  {   /* The previous line should work */
      char instStr[12];
      strcpy (instStr, "PESA HIGH");
      fprintf(stderr,"3D binning error: offset(%d) >= nsamples(%d)\n\r",
	      offset, nsamples);
      fprintf(stderr,"    ndata: %d, map type: %s\n\r", ndata, instStr);
    } else {
      if(offset < nsamples)
	alldata[offset++] = decomp19_8(*data++);
    }
  }
  status--;
  /* this is the **end** of the section executed on each call */
  
  if (status == 0) {   /* if all packets received, fill the data structure */

    uint ts=0;
    for(t=0;t<32;t++){ /* fill the phi-theta map */
      ts = tsect[t];
      for(p=0;p<32;p++){
        if(t<24)
        	map->pt_map[t*32+p] = ptmap[t*32+p];
	b = ptmap[ts*32+p];
	ps = (p - p0) & 0x1f;
	if (b > 87) 
	  printf("invalid map entry: t: %d, p: %d, ts: %d, ps: %d, t*32+ps: %d, b: %d\n\r",
		 t,p,ts,ps,t*32+ps,b);
      }
    }
    
    for(b=0;b<map->nbins;b++){
      offset         = bin[b].offset;
      map->domega[b] = bin[b].domega;
      for(e=0;e<map->nenergy;e++){
	map->nrg  [b][e] = nrg[e];
	map->dnrg [b][e] = dnrg[e];
	map->theta[b][e] = bin[b].theta;
        map->dtheta[b][e] = bin[b].dtheta;
        map->dphi  [b][e] = bin[b].dphi;
	map->gf   [b][e] = bin[b].geom;
	map->dt   [b][e] = num_phi[b] * dt;
	map->dvolume[b][e] = dnrg[e]/nrg[e]*map->domega[b];
	if (map->nenergy==bin[b].ne) {
	  flux  = alldata[offset+e];  
	  phi 	= bin[b].phi[e];
	} else if (map->nenergy==bin[b].ne/2) {
	  flux  = alldata[offset+2*e];
	  flux += alldata[offset+2*e+1];
	  phi	= bin[b].phi[2*e];
	  phi  += bin[b].phi[2*e+1];
	  phi	= phi/2.;
	} else {
	  flux  = 0;
	  /*printf("bin: %d: energy steps=%d, should be %d or %d\n\r",
		 b,bin[b].ne,map->nenergy,2*map->nenergy);*/
	}
	map->flux[b][e]  = flux;      }
 	map->phi[b][e] = phi;
   }

    if(double_sweep){
      int dsbins[4];
      int midbin;                       /* ramping bin between double sweeps */

      midbin=(map->nenergy-2+map->nenergy%2)/2;

      for(t=0;t<16;t++){
        dsbins[0] = ptmap[t*32+26];
        dsbins[1] = ptmap[t*32+27];
        dsbins[2] = ptmap[t*32+28];
        dsbins[3] = ptmap[t*32+29];
        if((dsbins[0] == dsbins[1]) || (dsbins[1] == dsbins[3]))
          for(i=0;i<6;i++)
	  for(e=midbin;e<map->nenergy;e++) {
	    map->flux[dsbins[1]][e] = NaN;        /* blank corrupted bins */
	    map->bins[dsbins[1]][e] = 0;        /* blank corrupted bins */
	  }
	else
	  for(i=0;i<16;i++){
	    map->nrg[dsbins[1]][midbin] = NaN;
	    map->dnrg[dsbins[1]][midbin] = NaN;
	    map->bins[dsbins[1]][midbin] = 0;        /* blank corrupted bins */
	    for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	      map->nrg [dsbins[1]][e] = map->nrg [dsbins[1]][e-midbin-1]; /*fix  e */
	      map->dnrg[dsbins[1]][e] = map->dnrg[dsbins[1]][e-midbin-1]; /*fix de */
	    }
	  }
        if(dsbins[1] != dsbins[2]) {
          if((dsbins[0] == dsbins[2]) || (dsbins[2] == dsbins[3]))
            for(i=0;i<6;i++)
  	    for(e=midbin;e<map->nenergy;e++) {
	      map->flux[dsbins[2]][e] = NaN;        /* blank corrupted bins */
	      map->bins[dsbins[2]][e] = 0;        /* blank corrupted bins */
	    }
	  else
	    for(i=0;i<16;i++){
	      map->nrg[dsbins[2]][midbin] = NaN;
	      map->dnrg[dsbins[2]][midbin] = NaN;
	      map->bins[dsbins[2]][midbin] = 0;        /* blank corrupted bins */
	      for(e=midbin+1;e<map->nenergy;e++){ /*e=7or14;e<15or30*/
	        map->nrg [dsbins[2]][e] = map->nrg [dsbins[2]][e-midbin-1]; /*fix  e */
	        map->dnrg[dsbins[2]][e] = map->dnrg[dsbins[2]][e-midbin-1]; /*fix de */
	      }
	    }
	}
      }
    }

    map->valid=1;
  }

  return(1);  
}

int ph15_to_idl(int argc, void *argv[]){
  int4 size,advance,index,*mapcode,*options,ok=0;
  uint4 ph_size;
  double *time;
  static packet_selector pks;
  idl_ph56  *idl56;
  idl_ph65  *idl65;
  idl_ph88  *idl88;
  idl_ph97  *idl97;  
  idl_ph121 *idl121; 

  options = (int4   *)argv[0];
  time    = (double *)argv[1];
  mapcode = (int4   *)argv[2];

  size    = options[0];
  advance = options[1];
  index   = options[2];

  /* negative index means get by time*/
  if      (advance) {SET_PKS_BY_INDEX(pks,pks.index+advance,P3D_ID);}
  else if (index<0) {SET_PKS_BY_TIME (pks,time[0],P3D_ID);          } 
  else              {SET_PKS_BY_INDEX(pks,index,P3D_ID);            }
  
  switch(*mapcode){
  case MAP11d: {
    ph_size = sizeof(idl_ph121);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl121 = (idl_ph121 *)argv[3];
    idl121->mapcode=*mapcode;
    ok = get_next_ph121_struct(&pks,idl121);
    idl121->index=pks.index;
    break; 
  }
  case MAP11b: {
    ph_size = sizeof(idl_ph97);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl97 = (idl_ph97  *)argv[3];
    idl97->mapcode=*mapcode;
    ok = get_next_ph97_struct (&pks,idl97);
    idl97->index=pks.index;
    break;
  }
  case MAP_8: {
    ph_size = sizeof(idl_ph56);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl56 = (idl_ph56  *)argv[3];
    idl56->mapcode=*mapcode;
    ok = get_next_ph56_struct (&pks,idl56);
    idl56->index=pks.index;
    break;
  }
  case MAP_0: {
    ph_size = sizeof(idl_ph65);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl65 = (idl_ph65  *)argv[3];
    idl65->mapcode=*mapcode;
    ok = get_next_ph65_struct (&pks,idl65);
    idl65->index=pks.index;
    break;
  }
  case MAP22d: {
    ph_size = sizeof(idl_ph88);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl88 = (idl_ph88  *)argv[3];
    idl88->mapcode=*mapcode;
    ok = get_next_ph88_struct (&pks,idl88);
    idl88->index=pks.index;
    break;
  }
  default:                                          
    printf("ph_dcm.c: Ph map type unknown: %d\n\r",*mapcode); /*ok=0*/
  }
  
  return(ok);
}

int ph30_to_idl(int argc, void *argv[]){ /* 30 energy steps, 5 angular bins */
  int4 size,advance,index,*mapcode,*options,ok=0;
  uint4 ph_size;
  double *time;
  static packet_selector pks;
  idl_ph5 *idl5;  

  options = (int4   *)argv[0];
  time    = (double *)argv[1];
  mapcode = (int4   *)argv[2];

  size    = options[0];
  advance = options[1];
  index   = options[2];

  /* negative index means get by time*/
  if      (advance) {SET_PKS_BY_INDEX(pks,pks.index+advance,P3D_ID);}
  else if (index<0) {SET_PKS_BY_TIME (pks,time[0],P3D_ID);          } 
  else              {SET_PKS_BY_INDEX(pks,index,P3D_ID);            }

  ph_size = sizeof(idl_ph5);
  if (size != ph_size) {
    printf("Incorrect input stucture size %d (should be %d).\r\n",
	   size,ph_size);
    return(0);
  }
  idl5 = (idl_ph5 *)argv[3];
  idl5->mapcode=*mapcode;
  ok = get_next_ph5_struct(&pks,idl5);  /* 5 angular bins, 30 energy steps */
  idl5->index=pks.index;
  
  return(ok);
}


int phb15_to_idl(int argc, void *argv[]){
  int4 size,advance,index,*mapcode,*options,ok=0;
  uint4 ph_size;
  double *time;
  static packet_selector pks;
  idl_ph56  *idl56;
  idl_ph65  *idl65;
  idl_ph88  *idl88;
  idl_ph97  *idl97;  
  idl_ph121 *idl121; 

  options = (int4   *)argv[0];
  time    = (double *)argv[1];
  mapcode = (int4   *)argv[2];

  size    = options[0];
  advance = options[1];
  index   = options[2];

  /* negative index means get by time*/
  if      (advance) {SET_PKS_BY_INDEX(pks,pks.index+advance,P3D_BRST_ID);}
  else if (index<0) {SET_PKS_BY_TIME (pks,time[0],P3D_BRST_ID);          } 
  else              {SET_PKS_BY_INDEX(pks,index,P3D_BRST_ID);            }

  switch(*mapcode){
  case MAP11d: {
    ph_size = sizeof(idl_ph121);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl121 = (idl_ph121 *)argv[3];
    idl121->mapcode=*mapcode;
    ok = get_next_ph121_struct(&pks,idl121);
    idl121->index=pks.index;
    break; 
  }
  case MAP11b: {
    ph_size = sizeof(idl_ph97);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl97 = (idl_ph97  *)argv[3];
    idl97->mapcode=*mapcode;
    ok = get_next_ph97_struct (&pks,idl97);
    idl97->index=pks.index;
    break;
  }
  case MAP_8: {
    ph_size = sizeof(idl_ph56);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl56 = (idl_ph56  *)argv[3];
    idl56->mapcode=*mapcode;
    ok = get_next_ph56_struct (&pks,idl56);
    idl56->index=pks.index;
    break;
  }
  case MAP_0: {
    ph_size = sizeof(idl_ph65);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl65 = (idl_ph65  *)argv[3];
    idl65->mapcode=*mapcode;
    ok = get_next_ph65_struct (&pks,idl65);
    idl65->index=pks.index;
    break;
  }
  case MAP22d: {
    ph_size = sizeof(idl_ph88);
    if (size != ph_size) {
      printf("Incorrect input stucture size %d (should be %d).\r\n",
	     size,ph_size);
      return(0);
    }
    idl88 = (idl_ph88  *)argv[3];
    idl88->mapcode=*mapcode;
    ok = get_next_ph88_struct (&pks,idl88);
    idl88->index=pks.index;
    break;
  }
  default:                                          
    printf("ph_dcm.c: Ph map type unknown: %d\n\r",*mapcode); /*ok=0*/
  }
  
  return(ok);
}

int4 ph15times_to_idl(int argc, void *argv[]){

  int4 size,advance,index,*mapcode,*options,ok=0;
  double *time;

  options = (int4   *)argv[0];
  time    = (double *)argv[1];

  size    = options[0];
  advance = options[1];
  index   = options[2];

  ok = get_time_points(P3D_ID,size,time);
  return(ok);
}

int4 phb15times_to_idl(int argc, void *argv[]){

  int4 size,advance,index,*mapcode,*options,ok=0;
  double *time;

  options = (int4   *)argv[0];
  time    = (double *)argv[1];

  size    = options[0];
  advance = options[1];
  index   = options[2];

  ok = get_time_points(P3D_BRST_ID,size,time);
  return(ok);
}