/*
 *  GVDR library for reading GVDR data files
 *  Copyright (C) 1994 Michael J. Maurer
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Library General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  This library 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
 *  Library General Public License for more details.
 *
 *  You should have received a copy of the GNU Library General Public
 *  License along with this library; if not, write to the Free
 *  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  Michael Maurer <maurer@nova.stanford.edu>
 *  Durand Bldg - Room 232
 *  Stanford, CA 94305-4055
 *  (415) 723-1024
 */
static char rcsid[]="$Id: tgvdrmat.c,v 1.7 1994/05/11 02:56:31 maurer Exp $";
/******************************************************************************
  tgvdrmat.c

     Function:  Matlab interface for tgvdr2mat.c.

This file is part of the STARLab Magellan Altimeter Data Processing Software.
Michael Maurer, August 1993.
******************************************************************************/
/* $Log: tgvdrmat.c,v $
 * Revision 1.7  1994/05/11  02:56:31  maurer
 * Changed calling syntax of tgvdrmat_init() and tgvdrmat_pixel().
 * Added match variable to speed things up when few variables are matched.
 *
 * Revision 1.6  1994/04/29  02:03:08  maurer
 * Updated to handle scattering law fits.
 *
 * Revision 1.5  1994/03/04  22:08:45  maurer
 * Now averages cohorts that match instead of taking first match only.
 *
 * Revision 1.4  1994/02/02  02:14:39  maurer
 * Changed ADF from min/max to mean/var.
 *
 * Revision 1.3  1994/01/06  00:21:50  maurer
 * Added ability to select cohorts based on certain fields.
 *
 * Revision 1.2  1993/11/02  21:17:25  maurer
 * Expanded to handle most GVDR data types.
 *
 * Revision 1.1  1993/09/10  19:43:23  maurer
 * Minor changes.
 *
 * Revision 1.0  1993/08/25  20:49:23  maurer
 * Initial revision
 * */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "libmisc.h"
#include "libml.h"
#include "gvdr.h"
#include "proj.h"	/* for gv_proj.p */
#include "tgvdrmat.h"

#define TGVDRMAT_C
#include "tgvdrmat.p"
#include "libml.p"
#include "gv_pack.p"	/* for cell_rlen() */
#include "gv_proj.p"

typedef unsigned char uchar;
typedef unsigned short ushort;

#ifdef HAVE_GNU_REGEX
static regex_t		s_pattern;
#else
int  re_exec();
char *re_comp();
#endif

#define DECLARE(type,var) \
   static type * var = NULL;\
   static float accum_##var;\
   static float lo_##var = -1e30;\
   static float hi_##var =  1e30;\
   static int match_##var

#define FITDECLARE(law) \
   DECLARE(uchar,gf_##law##_flags);\
   DECLARE(float,gf_##law##_p1);\
   DECLARE(float,gf_##law##_p1v);\
   DECLARE(float,gf_##law##_p2);\
   DECLARE(float,gf_##law##_p2v);\
   DECLARE(float,gf_##law##_rms);\
   DECLARE(float,gf_##law##_rmsv);\
   DECLARE(float,gf_##law##_resid);\

#ifdef HAVE_GNU_REGEX
#define MATCH(hit,var) do {\
   int i;\
   match_##var = !regexec(&s_pattern,#var,0,NULL,0);\
   if (match_##var) {\
      hit=1;\
      for (i=0; i<nre_lo; i++)\
	if (!regexec(&re_lo[i].pattern,#var,0,NULL,0)) {\
	   lo_##var = re_lo[i].val;\
	   break;\
	}\
      for (i=0; i<nre_hi; i++)\
	if (!regexec(&re_hi[i].pattern,#var,0,NULL,0)) {\
	   hi_##var = re_hi[i].val;\
	   break;\
	}\
   }\
} while (0)
#else
#define MATCH(hit,var) do {\
   int i;\
   match_##var = re_exec(#var);\
   if (match_##var) {\
      hit=1;\
      for (i=0; i<nre_lo; i++)\
	if (!strcmp(re_lo[i].regexp,#var)) {\
	   lo_##var = re_lo[i].val;\
	   break;\
	}\
      for (i=0; i<nre_hi; i++)\
	if (!strcmp(re_hi[i].regexp,#var)) {\
	   hi_##var = re_hi[i].val;\
	   break;\
	}\
   }\
} while (0)
#endif

#define FITMATCH(hit,law) do {\
   MATCH(hit,gf_##law##_flags);\
   MATCH(hit,gf_##law##_p1);\
   MATCH(hit,gf_##law##_p1v);\
   MATCH(hit,gf_##law##_p2);\
   MATCH(hit,gf_##law##_p2v);\
   MATCH(hit,gf_##law##_rms);\
   MATCH(hit,gf_##law##_rmsv);\
   MATCH(hit,gf_##law##_resid);\
} while (0)

#define ALLOC(type,var) do {\
   if (match_##var) \
      calloct_err((var),s_nx*s_ny,type);\
} while (0)

#define FITALLOC(law) do {\
   ALLOC(uchar,gf_##law##_flags);\
   ALLOC(float,gf_##law##_p1);\
   ALLOC(float,gf_##law##_p1v);\
   ALLOC(float,gf_##law##_p2);\
   ALLOC(float,gf_##law##_p2v);\
   ALLOC(float,gf_##law##_rms);\
   ALLOC(float,gf_##law##_rmsv);\
   ALLOC(float,gf_##law##_resid);\
} while (0)

#define ASSIGN0(var,xr,yr,src) do {\
   if (match_##var) \
      (var)[(yr)+s_ny*(xr)] = (src)->##var;\
} while (0)

#define FITASSIGN0(law,xr,yr,src) do {\
   ASSIGN0(gf_##law##_flags,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_p1,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_p1v,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_p2,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_p2v,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_rms,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_rmsv,	xr,yr,src,rec);\
   ASSIGN0(gf_##law##_resid,	xr,yr,src,rec);\
} while (0)

#define INRANGE(src,rec,coh,var) \
   ((src)->gv_##rec[coh].var >= lo_##var && (src)->gv_##rec[coh].var <= hi_##var)

#define ASSIGN(var,xr,yr,src,rec,coh,NaN) do {\
   if (match_##var)\
      (var)[(yr)+s_ny*(xr)] = (coh<(src)->gv_N##rec) ? \
         (src)->gv_##rec[coh].var : NaN;\
} while (0)

#define ASSIGNF(law,var,xr,yr,src,rec,coh,NaN) do {\
   if (match_gf_##law##_##var)\
      (var)[(yr)+s_ny*(xr)] = (coh<(src)->gv_N##rec && \
                               (src)->gv_##rec[coh].ga_fit[sl##law].gf_slaw == sl##law) ? \
         (src)->gv_##rec[coh].ga_fit[sl##law].gf_##var : NaN;\
} while (0)

#define FITASSIGN(law,xr,yr,src,rec,coh) do {\
   ASSIGNF(law,flags,	xr,yr,src,rec,coh,0);\
   ASSIGNF(law,p1,	xr,yr,src,rec,coh,s_NaN);\
   ASSIGNF(law,p1v,	xr,yr,src,rec,coh,s_NaN);\
   ASSIGNF(law,p2,	xr,yr,src,rec,coh,s_NaN);\
   ASSIGNF(law,p2v,	xr,yr,src,rec,coh,s_NaN);\
   ASSIGNF(law,rms,	xr,yr,src,rec,coh,s_NaN);\
   ASSIGNF(law,rmsv,	xr,yr,src,rec,coh,s_NaN);\
   ASSIGNF(law,resid,	xr,yr,src,rec,coh,s_NaN);\
} while (0)

#define NLOOKS0     s_Nlooks=0
#define ACCUM0(var) accum_##var=0.0

#define FITACCUM0(law) do {\
   ACCUM0(gf_##law##_flags);\
   ACCUM0(gf_##law##_p1);\
   ACCUM0(gf_##law##_p1v);\
   ACCUM0(gf_##law##_p2);\
   ACCUM0(gf_##law##_p2v);\
   ACCUM0(gf_##law##_rms);\
   ACCUM0(gf_##law##_rmsv);\
   ACCUM0(gf_##law##_resid);\
} while (0)

#define NLOOKSINC(var,xr,yr,src,rec,coh) do {\
   s_Nlooks += (src)->gv_##rec[coh].var;\
   if (match_##var)\
      (var)[(yr)+s_ny*(xr)] = s_Nlooks;\
} while (0)

#define ASSIGNINC(var,xr,yr,src,rec,coh,inc) do {\
   if (match_##var)\
      accum_##var += (src)->gv_##rec[coh].inc*(src)->gv_##rec[coh].var;\
} while (0)

#define ASSIGNINCF(law,var,xr,yr,src,rec,coh,inc) do {\
   if (match_gf_##law##_##var)\
      accum_gf_##law##_##var += (src)->gv_##rec[coh].inc * \
                                (src)->gv_##rec[coh].ga_fit[sl##law].gf_##var;\
} while (0)

#define FITASSIGNINC(law,xr,yr,src,rec,coh,inc) do {\
   ASSIGNINCF(law,flags,xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,p1,	xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,p1v,	xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,p2,	xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,p2v,	xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,rms,	xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,rmsv,	xr,yr,src,rec,coh,inc);\
   ASSIGNINCF(law,resid,xr,yr,src,rec,coh,inc);\
} while (0)

#define ASSIGNDIV(var,xr,yr,NaN) do {\
   if (match_##var)\
      (var)[(yr)+s_ny*(xr)] = s_Nlooks ? \
                accum_##var / s_Nlooks : \
                NaN;\
} while (0)

#define FITASSIGNDIV(law,xr,yr) do {\
   ASSIGNDIV(gf_##law##_flags,	xr,yr,0);\
   ASSIGNDIV(gf_##law##_p1,	xr,yr,s_NaN);\
   ASSIGNDIV(gf_##law##_p1v,	xr,yr,s_NaN);\
   ASSIGNDIV(gf_##law##_p2,	xr,yr,s_NaN);\
   ASSIGNDIV(gf_##law##_p2v,	xr,yr,s_NaN);\
   ASSIGNDIV(gf_##law##_rms,	xr,yr,s_NaN);\
   ASSIGNDIV(gf_##law##_rmsv,	xr,yr,s_NaN);\
   ASSIGNDIV(gf_##law##_resid,	xr,yr,s_NaN);\
} while (0)

#define WRITE(fout,var,mtype) do {\
   if (match_##var) {\
      mlwrite(fout,ML_TYPE+(mtype),#var,s_ny,s_nx,0,(var),NULL);\
      free(var);\
      (var)=NULL;\
   }\
} while (0)

#define FITWRITE(fout,law) do {\
   WRITE(fout,gf_##law##_flags,	ML_UCHAR);\
   WRITE(fout,gf_##law##_p1,	ML_FLOAT);\
   WRITE(fout,gf_##law##_p1v,	ML_FLOAT);\
   WRITE(fout,gf_##law##_p2,	ML_FLOAT);\
   WRITE(fout,gf_##law##_p2v,	ML_FLOAT);\
   WRITE(fout,gf_##law##_rms,	ML_FLOAT);\
   WRITE(fout,gf_##law##_rmsv,	ML_FLOAT);\
   WRITE(fout,gf_##law##_resid,	ML_FLOAT);\
} while (0)

static int	s_init=0;
static double	s_NaN;
static int	s_nx;
static int	s_ny;
static proj_t	s_P;
static int	s_Nlooks;

DECLARE(int,	gv_X);
DECLARE(int,	gv_Y);
DECLARE(double,	gv_lon);
DECLARE(double,	gv_lat);
DECLARE(int,	gv_reclen);
DECLARE(uchar,	gv_Nxif);
DECLARE(uchar,	gv_Nedf);
DECLARE(uchar,	gv_Nadf);
DECLARE(uchar,	gv_Nanf);

DECLARE(ushort,	gx_Nlooks);
DECLARE(float,	gx_azim);
DECLARE(float,	gx_inc);
DECLARE(float,	gx_poln);
DECLARE(uchar,	gx_hlo);
DECLARE(uchar,	gx_hmed);
DECLARE(uchar,	gx_hhi);
DECLARE(uchar,	gx_hpk);
DECLARE(float,	gx_c0);
DECLARE(float,	gx_c1);
DECLARE(float,	gx_c2);

DECLARE(ushort,	ge_Nlooks);
DECLARE(float,	ge_azim);
DECLARE(float,	ge_inc);
DECLARE(float,	ge_poln);
DECLARE(float,	ge_emiss);
DECLARE(float,	ge_emissv);

DECLARE(ushort,	gd_Nlooks);
DECLARE(float,	gd_radius);
DECLARE(float,	gd_radiusv);
DECLARE(float,	gd_slope);
DECLARE(float,	gd_slopev);
DECLARE(float,	gd_rho);
DECLARE(float,	gd_rhov);

DECLARE(ushort,	ga_Nlooks);
DECLARE(uchar,	ga_Nsig0);
DECLARE(uchar,	ga_Nfit);
DECLARE(float,	ga_dcent);
DECLARE(float,	ga_azim);
/* skip sig0, Vsig0 */

FITDECLARE(Hag);
FITDECLARE(Exp);
FITDECLARE(Gau);
FITDECLARE(Ray);
FITDECLARE(Muh);

int
tgvdrmat_init(gvdr_sfduhdr_t *gh,
	      int nx,
	      int ny,
	      char *regexp,
	      re_limit_t *re_lo,
	      re_limit_t *re_hi,
	      int nre_lo,
	      int nre_hi,
	      gvmatch_t *match)
{
   int err,i;
   char msg1[256],*msg2;
   double sqrt();

   s_nx=nx;
   s_ny=ny;
#ifdef HAVE_GNU_REGEX
   /* use GNU POSIX matching */
   if (err=regcomp(&s_pattern,regexp,REG_EXTENDED)) {
      regerror(err,&s_pattern,msg1,sizeof(msg1));
      fprintf(stderr,"%s: regcomp failed on %s: %s\n",progn,regexp,msg1);
      return 1;
   }
   for (i=0; i<nre_lo; i++)
      if (err=regcomp(&re_lo[i].pattern,re_lo[i].regexp,REG_EXTENDED)) {
      regerror(err,&re_lo[i].pattern,msg1,sizeof(msg1));
      fprintf(stderr,"%s: regcomp failed on %s: %s\n",progn,re_lo[i].regexp,msg1);
      return 2;
   }
   for (i=0; i<nre_hi; i++)
      if (err=regcomp(&re_hi[i].pattern,re_hi[i].regexp,REG_EXTENDED)) {
      regerror(err,&re_hi[i].pattern,msg1,sizeof(msg1));
      fprintf(stderr,"%s: regcomp failed on %s: %s\n",progn,re_hi[i].regexp,msg1);
      return 2;
   }
#else
   /* use ed(1) style matching */
   if (msg2=re_comp(regexp)) {
      fprintf(stderr,"%s: re_comp failed on %s: %s\n",progn,regexp,msg2);
      return 1;
   }
#endif

   match->gv = match->gx = match->ge = match->gd = match->ga = 0;

   MATCH(match->xy,gv_X);
   MATCH(match->xy,gv_Y);
   MATCH(match->xy,gv_lon);
   MATCH(match->xy,gv_lat);
   MATCH(match->gv,gv_reclen);
   MATCH(match->gv,gv_Nxif);
   MATCH(match->gv,gv_Nedf);
   MATCH(match->gv,gv_Nadf);
   MATCH(match->gv,gv_Nanf);

   MATCH(match->gx,gx_Nlooks);
   MATCH(match->gx,gx_azim);
   MATCH(match->gx,gx_inc);
   MATCH(match->gx,gx_poln);
   MATCH(match->gx,gx_hlo);
   MATCH(match->gx,gx_hmed);
   MATCH(match->gx,gx_hhi);
   MATCH(match->gx,gx_hpk);
   MATCH(match->gx,gx_c0);
   MATCH(match->gx,gx_c1);
   MATCH(match->gx,gx_c2);

   MATCH(match->ge,ge_Nlooks);
   MATCH(match->ge,ge_azim);
   MATCH(match->ge,ge_inc);
   MATCH(match->ge,ge_poln);
   MATCH(match->ge,ge_emiss);
   MATCH(match->ge,ge_emissv);

   MATCH(match->gd,gd_Nlooks);
   MATCH(match->gd,gd_radius);
   MATCH(match->gd,gd_radiusv);
   MATCH(match->gd,gd_slope);
   MATCH(match->gd,gd_slopev);
   MATCH(match->gd,gd_rho);
   MATCH(match->gd,gd_rhov);

   MATCH(match->ga,ga_Nlooks);
   MATCH(match->ga,ga_Nsig0);
   MATCH(match->ga,ga_Nfit);
   MATCH(match->ga,ga_dcent);
   MATCH(match->ga,ga_azim);

   FITMATCH(match->ga,Hag);
   FITMATCH(match->ga,Exp);
   FITMATCH(match->ga,Gau);
   FITMATCH(match->ga,Ray);
   FITMATCH(match->ga,Muh);

   match->dat = match->gv || match->gx || match->ge || match->gd || match->ga;
   match->any = match->xy || match->dat;

   ALLOC(int,	gv_X);
   ALLOC(int,	gv_Y);
   ALLOC(double,gv_lon);
   ALLOC(double,gv_lat);
   ALLOC(int,	gv_reclen);
   ALLOC(uchar,	gv_Nxif);
   ALLOC(uchar,	gv_Nedf);
   ALLOC(uchar,	gv_Nadf);
   ALLOC(uchar,	gv_Nanf);

   ALLOC(ushort,gx_Nlooks);
   ALLOC(float,	gx_azim);
   ALLOC(float,	gx_inc);
   ALLOC(float,	gx_poln);
   ALLOC(uchar,	gx_hlo);
   ALLOC(uchar,	gx_hmed);
   ALLOC(uchar,	gx_hhi);
   ALLOC(uchar,	gx_hpk);
   ALLOC(float,	gx_c0);
   ALLOC(float,	gx_c1);
   ALLOC(float,	gx_c2);

   ALLOC(ushort,ge_Nlooks);
   ALLOC(float,	ge_azim);
   ALLOC(float,	ge_inc);
   ALLOC(float,	ge_poln);
   ALLOC(float,	ge_emiss);
   ALLOC(float,	ge_emissv);

   ALLOC(ushort,gd_Nlooks);
   ALLOC(float,	gd_radius);
   ALLOC(float,	gd_radiusv);
   ALLOC(float,	gd_slope);
   ALLOC(float,	gd_slopev);
   ALLOC(float,	gd_rho);
   ALLOC(float,	gd_rhov);

   ALLOC(ushort,ga_Nlooks);
   ALLOC(uchar,	ga_Nsig0);
   ALLOC(uchar,	ga_Nfit);
   ALLOC(float,	ga_dcent);
   ALLOC(float,	ga_azim);

   FITALLOC(Hag);
   FITALLOC(Exp);
   FITALLOC(Gau);
   FITALLOC(Ray);
   FITALLOC(Muh);

   s_NaN=sqrt(-1.0);

   if (match_gv_lat || match_gv_lon) {
      if (err=gvdr_rproj(gh,&s_P)) {
	 error(0,0,"[tgvdrmat_init] gvdr_rproj failed %d",err);
	 return 1;
      }
   }

   s_init=1;
   return 0;
}

int
tgvdrmat_pixel(int xr,
	       int yr,
	       gvdrcell_t *C,
	       gvmatch_t *match)
{
   if (!s_init)
      error(-1,0,"[tgvdrmat_pixel] package not initialized");

   if (xr<0 || xr>=s_nx || yr<0 || yr>=s_ny)
      error(-1,0,"[tgvdrmat_pixel] xr=%d yr=%d out of range",xr,yr);

   /* fill otherwise empty fields if needed */
   C->gv_reclen=cell_rlen(C);
   if (match_gv_lat || match_gv_lon)
      gvdr_tolat(&s_P,C);

   if (match->xy || match->gv) tgvdrmat_hdr(xr,yr,C);
   if (match->gx) tgvdrmat_xif(xr,yr,C);
   if (match->ge) tgvdrmat_edf(xr,yr,C);
   if (match->gd) tgvdrmat_adf(xr,yr,C);
   if (match->ga) tgvdrmat_anf(xr,yr,C);
   return 0;
}

int
tgvdrmat_hdr(int xr,
	     int yr,
	     gvdrcell_t *C)
{
   ASSIGN0(gv_X,	xr,yr,C);
   ASSIGN0(gv_Y,	xr,yr,C);
   ASSIGN0(gv_lon,	xr,yr,C);
   ASSIGN0(gv_lat,	xr,yr,C);
   ASSIGN0(gv_reclen,	xr,yr,C);
   ASSIGN0(gv_Nxif,	xr,yr,C);
   ASSIGN0(gv_Nedf,	xr,yr,C);
   ASSIGN0(gv_Nadf,	xr,yr,C);
   ASSIGN0(gv_Nanf,	xr,yr,C);
   return 0;
}

int
tgvdrmat_xif(int xr,
	     int yr,
	     gvdrcell_t *C)
{
   int coh;

   NLOOKS0;
   ACCUM0(gx_azim);
   ACCUM0(gx_inc);
   ACCUM0(gx_poln);
   ACCUM0(gx_hlo);
   ACCUM0(gx_hmed);
   ACCUM0(gx_hhi);
   ACCUM0(gx_hpk);
   ACCUM0(gx_c0);
   ACCUM0(gx_c1);
   ACCUM0(gx_c2);
   for (coh=0; coh<C->gv_Nxif; coh++) {
      if (INRANGE(C,xif,coh,gx_Nlooks) &&
	  INRANGE(C,xif,coh,gx_azim) &&
	  INRANGE(C,xif,coh,gx_inc) &&
	  INRANGE(C,xif,coh,gx_poln) &&
	  INRANGE(C,xif,coh,gx_hlo) &&
	  INRANGE(C,xif,coh,gx_hmed) &&
	  INRANGE(C,xif,coh,gx_hhi) &&
	  INRANGE(C,xif,coh,gx_hpk) &&
	  INRANGE(C,xif,coh,gx_c0) &&
	  INRANGE(C,xif,coh,gx_c1) &&
	  INRANGE(C,xif,coh,gx_c2)) {
	 NLOOKSINC(gx_Nlooks,	xr,yr,C,xif,coh);
	 ASSIGNINC(gx_azim,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_inc,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_poln,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_hlo,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_hmed,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_hhi,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_hpk,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_c0,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_c1,	xr,yr,C,xif,coh,gx_Nlooks);
	 ASSIGNINC(gx_c2,	xr,yr,C,xif,coh,gx_Nlooks);
      }
   }
   ASSIGNDIV(gx_azim,	xr,yr,s_NaN);
   ASSIGNDIV(gx_inc,	xr,yr,s_NaN);
   ASSIGNDIV(gx_poln,	xr,yr,s_NaN);
   ASSIGNDIV(gx_hlo,	xr,yr,255);
   ASSIGNDIV(gx_hmed,	xr,yr,255);
   ASSIGNDIV(gx_hhi,	xr,yr,255);
   ASSIGNDIV(gx_hpk,	xr,yr,255);
   ASSIGNDIV(gx_c0,	xr,yr,s_NaN);
   ASSIGNDIV(gx_c1,	xr,yr,s_NaN);
   ASSIGNDIV(gx_c2,	xr,yr,s_NaN);
   return 0;
}

int
tgvdrmat_edf(int xr,
	     int yr,
	     gvdrcell_t *C)
{
   int coh;

   NLOOKS0;
   ACCUM0(ge_azim);
   ACCUM0(ge_inc);
   ACCUM0(ge_poln);
   ACCUM0(ge_emiss);
   ACCUM0(ge_emissv);
   for (coh=0; coh<C->gv_Nedf; coh++) {
      if (INRANGE(C,edf,coh,ge_Nlooks) &&
	  INRANGE(C,edf,coh,ge_azim) &&
	  INRANGE(C,edf,coh,ge_inc) &&
	  INRANGE(C,edf,coh,ge_poln) &&
	  INRANGE(C,edf,coh,ge_emiss) &&
	  INRANGE(C,edf,coh,ge_emissv)) {
	 NLOOKSINC(ge_Nlooks,	xr,yr,C,edf,coh);
	 ASSIGNINC(ge_azim,	xr,yr,C,edf,coh,ge_Nlooks);
	 ASSIGNINC(ge_inc,	xr,yr,C,edf,coh,ge_Nlooks);
	 ASSIGNINC(ge_poln,	xr,yr,C,edf,coh,ge_Nlooks);
	 ASSIGNINC(ge_emiss,	xr,yr,C,edf,coh,ge_Nlooks);
	 ASSIGNINC(ge_emissv,	xr,yr,C,edf,coh,ge_Nlooks);
      }
   }
   ASSIGNDIV(ge_azim,	xr,yr,s_NaN);
   ASSIGNDIV(ge_inc,	xr,yr,s_NaN);
   ASSIGNDIV(ge_poln,	xr,yr,s_NaN);
   ASSIGNDIV(ge_emiss,	xr,yr,s_NaN);
   ASSIGNDIV(ge_emissv,	xr,yr,s_NaN);
   return 0;
}

int
tgvdrmat_adf(int xr,
	     int yr,
	     gvdrcell_t *C)
{
   int coh;

   NLOOKS0;
   ACCUM0(gd_radius);
   ACCUM0(gd_radiusv);
   ACCUM0(gd_slope);
   ACCUM0(gd_slopev);
   ACCUM0(gd_rho);
   ACCUM0(gd_rhov);
   for (coh=0; coh<C->gv_Nadf; coh++) {
      if (INRANGE(C,adf,coh,gd_Nlooks) &&
	  INRANGE(C,adf,coh,gd_radius) &&
	  INRANGE(C,adf,coh,gd_radiusv) &&
	  INRANGE(C,adf,coh,gd_slope) &&
	  INRANGE(C,adf,coh,gd_slopev) &&
	  INRANGE(C,adf,coh,gd_rho) &&
	  INRANGE(C,adf,coh,gd_rhov)) {
	 NLOOKSINC(gd_Nlooks,	xr,yr,C,adf,coh);
	 ASSIGNINC(gd_radius,	xr,yr,C,adf,coh,gd_Nlooks);
	 ASSIGNINC(gd_radiusv,	xr,yr,C,adf,coh,gd_Nlooks);
	 ASSIGNINC(gd_slope,	xr,yr,C,adf,coh,gd_Nlooks);
	 ASSIGNINC(gd_slopev,	xr,yr,C,adf,coh,gd_Nlooks);
	 ASSIGNINC(gd_rho,	xr,yr,C,adf,coh,gd_Nlooks);
	 ASSIGNINC(gd_rhov,	xr,yr,C,adf,coh,gd_Nlooks);
      }
   }
   ASSIGNDIV(gd_radius,	xr,yr,s_NaN);
   ASSIGNDIV(gd_radiusv,xr,yr,s_NaN);
   ASSIGNDIV(gd_slope,	xr,yr,s_NaN);
   ASSIGNDIV(gd_slopev,	xr,yr,s_NaN);
   ASSIGNDIV(gd_rho,	xr,yr,s_NaN);
   ASSIGNDIV(gd_rhov,	xr,yr,s_NaN);
   return 0;
}

int
tgvdrmat_anf(int xr,
	     int yr,
	     gvdrcell_t *C)
{
   int coh;

   NLOOKS0;
   ACCUM0(ga_Nsig0);
   ACCUM0(ga_Nfit);
   ACCUM0(ga_dcent);
   ACCUM0(ga_azim);
   FITACCUM0(Hag);
   FITACCUM0(Exp);
   FITACCUM0(Gau);
   FITACCUM0(Ray);
   FITACCUM0(Muh);
   for (coh=0; coh<C->gv_Nanf; coh++) {
      if (INRANGE(C,anf,coh,ga_Nlooks) &&
	  INRANGE(C,anf,coh,ga_Nsig0) &&
	  INRANGE(C,anf,coh,ga_Nfit) &&
	  INRANGE(C,anf,coh,ga_dcent) &&
	  INRANGE(C,anf,coh,ga_azim)) {
	 NLOOKSINC(ga_Nlooks,	xr,yr,C,anf,coh);
	 ASSIGNINC(ga_Nsig0,	xr,yr,C,anf,coh,ga_Nlooks);
	 ASSIGNINC(ga_Nfit,	xr,yr,C,anf,coh,ga_Nlooks);
	 ASSIGNINC(ga_dcent,	xr,yr,C,anf,coh,ga_Nlooks);
	 ASSIGNINC(ga_azim,	xr,yr,C,anf,coh,ga_Nlooks);
	 FITASSIGNINC(Hag,	xr,yr,C,anf,coh,ga_Nlooks);
	 FITASSIGNINC(Exp,	xr,yr,C,anf,coh,ga_Nlooks);
	 FITASSIGNINC(Gau,	xr,yr,C,anf,coh,ga_Nlooks);
	 FITASSIGNINC(Ray,	xr,yr,C,anf,coh,ga_Nlooks);
	 FITASSIGNINC(Muh,	xr,yr,C,anf,coh,ga_Nlooks);
      }
   }
   ASSIGNDIV(ga_Nsig0,	xr,yr,0);
   ASSIGNDIV(ga_Nfit,	xr,yr,0);
   ASSIGNDIV(ga_dcent,	xr,yr,s_NaN);
   ASSIGNDIV(ga_azim,	xr,yr,s_NaN);
   FITASSIGNDIV(Hag,	xr,yr);
   FITASSIGNDIV(Exp,	xr,yr);
   FITASSIGNDIV(Gau,	xr,yr);
   FITASSIGNDIV(Ray,	xr,yr);
   FITASSIGNDIV(Muh,	xr,yr);
   return 0;
}

int
tgvdrmat_finish(FILE *fout)
{

   WRITE(fout,gv_X,		ML_INT);
   WRITE(fout,gv_Y,		ML_INT);
   WRITE(fout,gv_lon,		ML_DOUBLE);
   WRITE(fout,gv_lat,		ML_DOUBLE);
   WRITE(fout,gv_reclen,	ML_INT);
   WRITE(fout,gv_Nxif,		ML_UCHAR);
   WRITE(fout,gv_Nedf,		ML_UCHAR);
   WRITE(fout,gv_Nadf,		ML_UCHAR);
   WRITE(fout,gv_Nanf,		ML_UCHAR);

   WRITE(fout,gx_Nlooks,	ML_USHORT);
   WRITE(fout,gx_azim,		ML_FLOAT);
   WRITE(fout,gx_inc,		ML_FLOAT);
   WRITE(fout,gx_poln,		ML_FLOAT);
   WRITE(fout,gx_hlo,		ML_UCHAR);
   WRITE(fout,gx_hmed,		ML_UCHAR);
   WRITE(fout,gx_hhi,		ML_UCHAR);
   WRITE(fout,gx_hpk,		ML_UCHAR);
   WRITE(fout,gx_c0,		ML_FLOAT);
   WRITE(fout,gx_c1,		ML_FLOAT);
   WRITE(fout,gx_c2,		ML_FLOAT);

   WRITE(fout,ge_Nlooks,	ML_USHORT);
   WRITE(fout,ge_azim,		ML_FLOAT);
   WRITE(fout,ge_inc,		ML_FLOAT);
   WRITE(fout,ge_poln,		ML_FLOAT);
   WRITE(fout,ge_emiss,		ML_FLOAT);
   WRITE(fout,ge_emissv,	ML_FLOAT);

   WRITE(fout,gd_Nlooks,	ML_USHORT);
   WRITE(fout,gd_radius,	ML_FLOAT);
   WRITE(fout,gd_radiusv,	ML_FLOAT);
   WRITE(fout,gd_slope,		ML_FLOAT);
   WRITE(fout,gd_slopev,	ML_FLOAT);
   WRITE(fout,gd_rho,		ML_FLOAT);
   WRITE(fout,gd_rhov,		ML_FLOAT);

   WRITE(fout,ga_Nlooks,	ML_USHORT);
   WRITE(fout,ga_Nsig0,		ML_UCHAR);
   WRITE(fout,ga_Nfit,		ML_UCHAR);
   WRITE(fout,ga_dcent,		ML_FLOAT);
   WRITE(fout,ga_azim,		ML_FLOAT);

   FITWRITE(fout,Hag);
   FITWRITE(fout,Exp);
   FITWRITE(fout,Gau);
   FITWRITE(fout,Ray);
   FITWRITE(fout,Muh);

   s_init=0;
   return 0;
}