core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Add assertion and clean up indents

This commit is contained in:
Shahram Najm 2013-07-23 17:21:08 +01:00
parent 3d6ebcfaf3
commit 6c03d54a38
9 changed files with 1032 additions and 1030 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

152
src/grib_accessor_class_codeflag.c
View File

@ -133,108 +133,108 @@ static void init_class(grib_accessor_class* c)
static void init (grib_accessor* a,const long len, grib_arguments* param) static void init (grib_accessor* a,const long len, grib_arguments* param)
{ {
grib_accessor_codeflag* self = (grib_accessor_codeflag*)a; grib_accessor_codeflag* self = (grib_accessor_codeflag*)a;
a->length = len; a->length = len;
self->tablename = grib_arguments_get_string(a->parent->h,param,0); self->tablename = grib_arguments_get_string(a->parent->h,param,0);
Assert(a->length>=0); Assert(a->length>=0);
} }
static int test_bit(long a, long b) static int test_bit(long a, long b)
{ {
return a&(1<<b); return a&(1<<b);
} }
static int grib_get_codeflag(grib_accessor* a, long code, char* codename) static int grib_get_codeflag(grib_accessor* a, long code, char* codename)
{ {
grib_accessor_codeflag* self = (grib_accessor_codeflag*)a; grib_accessor_codeflag* self = (grib_accessor_codeflag*)a;
FILE* f = NULL; FILE* f = NULL;
char fname[1024]; char fname[1024];
char bval[50]; char bval[50];
char num[50]; char num[50];
char* filename=0; char* filename=0;
char line[1024]; char line[1024];
size_t i =0; size_t i =0;
int j =0; int j =0;
grib_recompose_name(a->parent->h,NULL, self->tablename, fname,1); grib_recompose_name(a->parent->h,NULL, self->tablename, fname,1);
if ((filename=grib_context_full_defs_path(a->parent->h->context,fname))==NULL) { if ((filename=grib_context_full_defs_path(a->parent->h->context,fname))==NULL) {
grib_context_log(a->parent->h->context,GRIB_LOG_WARNING,"Cannot open flag table %s",filename); grib_context_log(a->parent->h->context,GRIB_LOG_WARNING,"Cannot open flag table %s",filename);
strcpy(codename, "Cannot open flag table"); strcpy(codename, "Cannot open flag table");
return GRIB_FILE_NOT_FOUND; return GRIB_FILE_NOT_FOUND;
}
f=fopen(filename, "r");
if (!f)
{
grib_context_log(a->parent->h->context,(GRIB_LOG_WARNING)|(GRIB_LOG_PERROR),"Cannot open flag table %s",filename);
strcpy(codename, "Cannot open flag table");
return GRIB_FILE_NOT_FOUND;
}
#if 0
strcpy(codename, self->tablename);
strcat(codename,": ");
j = strlen(codename);
#endif
while(fgets(line,sizeof(line)-1,f))
{
sscanf(line,"%s %s", num, bval);
if(num[0] != '#')
{
if((test_bit(code, a->length*8-atol(num))>0) == atol(bval))
{
size_t linelen = strlen(line);
codename[j++] = '(';
codename[j++] = num[0];
codename[j++] = '=';
codename[j++] = bval[0];
codename[j++] = ')';
codename[j++] = ' ';
if(j)
codename[j++] = ' ';
for(i=(strlen(num)+strlen(bval)+2); i < linelen-1;i++)
codename[j++] = line[i];
if(line[i]!='\n')
codename[j++] = line[i];
codename[j++] = ';';
}
}
} }
if(j>1 && codename[j-1] == ';') j--; f=fopen(filename, "r");
codename[j] = 0;
strcat(codename,":"); if (!f)
strcat(codename,self->tablename); {
grib_context_log(a->parent->h->context,(GRIB_LOG_WARNING)|(GRIB_LOG_PERROR),"Cannot open flag table %s",filename);
strcpy(codename, "Cannot open flag table");
return GRIB_FILE_NOT_FOUND;
}
fclose(f); #if 0
return GRIB_SUCCESS; strcpy(codename, self->tablename);
strcat(codename,": ");
j = strlen(codename);
#endif
while(fgets(line,sizeof(line)-1,f))
{
sscanf(line,"%s %s", num, bval);
if(num[0] != '#')
{
if((test_bit(code, a->length*8-atol(num))>0) == atol(bval))
{
size_t linelen = strlen(line);
codename[j++] = '(';
codename[j++] = num[0];
codename[j++] = '=';
codename[j++] = bval[0];
codename[j++] = ')';
codename[j++] = ' ';
if(j)
codename[j++] = ' ';
for(i=(strlen(num)+strlen(bval)+2); i < linelen-1;i++)
codename[j++] = line[i];
if(line[i]!='\n')
codename[j++] = line[i];
codename[j++] = ';';
}
}
}
if(j>1 && codename[j-1] == ';') j--;
codename[j] = 0;
strcat(codename,":");
strcat(codename,self->tablename);
fclose(f);
return GRIB_SUCCESS;
} }
static long value_count(grib_accessor* a) static long value_count(grib_accessor* a)
{ {
return 1; return 1;
} }
static void dump(grib_accessor* a, grib_dumper* dumper) static void dump(grib_accessor* a, grib_dumper* dumper)
{ {
grib_accessor_codeflag* self = (grib_accessor_codeflag*)a; grib_accessor_codeflag* self = (grib_accessor_codeflag*)a;
long v; long v;
char flagname[1024]; char flagname[1024];
char fname[1024]; char fname[1024];
size_t llen = 1; size_t llen = 1;
grib_recompose_name(a->parent->h,NULL, self->tablename, fname,1); grib_recompose_name(a->parent->h,NULL, self->tablename, fname,1);
grib_unpack_long(a, &v, &llen); grib_unpack_long(a, &v, &llen);
grib_get_codeflag(a, v, flagname); grib_get_codeflag(a, v, flagname);
grib_dump_bits(dumper,a,flagname); grib_dump_bits(dumper,a,flagname);
} }

124
src/grib_accessor_class_g1bitmap.c
View File

@ -136,100 +136,100 @@ static void grib_set_bit_on( unsigned char* p, long *bitp);
static void init(grib_accessor* a, const long len , grib_arguments* arg ) static void init(grib_accessor* a, const long len , grib_arguments* arg )
{ {
grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a; grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a;
self->unusedBits = grib_arguments_get_name(a->parent->h,arg,4); self->unusedBits = grib_arguments_get_name(a->parent->h,arg,4);
} }
static int pack_double(grib_accessor* a, const double* val,size_t *len){ static int pack_double(grib_accessor* a, const double* val,size_t *len){
grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a; grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a;
size_t tlen; size_t tlen;
unsigned char* buf = NULL; unsigned char* buf = NULL;
size_t i; size_t i;
int err = 0; int err = 0;
long pos = 0; long pos = 0;
long bmaplen = 0; long bmaplen = 0;
const int bit_padding = 16; const int bit_padding = 16;
double miss_values = 0; double miss_values = 0;
tlen = ((*len+bit_padding-1)/bit_padding*bit_padding)/8; tlen = ((*len+bit_padding-1)/bit_padding*bit_padding)/8;
if((err = grib_get_double_internal(a->parent->h, self->missing_value, &miss_values)) if((err = grib_get_double_internal(a->parent->h, self->missing_value, &miss_values))
!= GRIB_SUCCESS) != GRIB_SUCCESS)
return err; return err;
buf = grib_context_malloc_clear(a->parent->h->context,tlen); buf = grib_context_malloc_clear(a->parent->h->context,tlen);
if(!buf) return GRIB_OUT_OF_MEMORY; if(!buf) return GRIB_OUT_OF_MEMORY;
pos=0; pos=0;
for(i=0;i<*len;i++) for(i=0;i<*len;i++)
{ {
if (val[i] == miss_values) if (val[i] == miss_values)
pos++; pos++;
else{ else{
bmaplen++; bmaplen++;
grib_set_bit_on(buf, &pos); grib_set_bit_on(buf, &pos);
}
} }
}
if((err = grib_set_long_internal(a->parent->h, self->unusedBits,tlen*8 - *len )) if((err = grib_set_long_internal(a->parent->h, self->unusedBits,tlen*8 - *len ))
!= GRIB_SUCCESS) != GRIB_SUCCESS)
return err; return err;
grib_buffer_replace(a, buf, tlen,1,1); grib_buffer_replace(a, buf, tlen,1,1);
grib_context_free(a->parent->h->context,buf); grib_context_free(a->parent->h->context,buf);
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static long value_count(grib_accessor* a) static long value_count(grib_accessor* a)
{ {
grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a; grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a;
long tlen; long tlen;
int err; int err;
if ((err=grib_get_long_internal(a->parent->h, self->unusedBits, &tlen)) != GRIB_SUCCESS) if ((err=grib_get_long_internal(a->parent->h, self->unusedBits, &tlen)) != GRIB_SUCCESS)
grib_context_log(a->parent->h->context, GRIB_LOG_ERROR, "grib_accessor_class_bitmap.value_count : cannot get %s err=%d",self->unusedBits,err); grib_context_log(a->parent->h->context, GRIB_LOG_ERROR, "grib_accessor_class_bitmap.value_count : cannot get %s err=%d",self->unusedBits,err);
tlen = (a->length*8)-tlen; tlen = (a->length*8)-tlen;
return tlen; return tlen;
} }
static int unpack_bytes(grib_accessor* a, unsigned char* val, size_t *len) static int unpack_bytes(grib_accessor* a, unsigned char* val, size_t *len)
{ {
unsigned char* buf = a->parent->h->buffer->data; unsigned char* buf = a->parent->h->buffer->data;
grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a; grib_accessor_g1bitmap* self = (grib_accessor_g1bitmap*)a;
long tlen; long tlen;
int err; int err;
long length = grib_byte_count(a); long length = grib_byte_count(a);
long offset = grib_byte_offset(a); long offset = grib_byte_offset(a);
if(*len < (size_t)length ) if(*len < (size_t)length )
{ {
grib_context_log(a->parent->h->context, GRIB_LOG_ERROR, "Wrong size for %s it is %d bytes long\n", a->name ,length ); grib_context_log(a->parent->h->context, GRIB_LOG_ERROR, "Wrong size for %s it is %d bytes long\n", a->name ,length );
*len = length;
return GRIB_ARRAY_TOO_SMALL;
}
if ((err=grib_get_long_internal(a->parent->h, self->unusedBits, &tlen)) != GRIB_SUCCESS)
grib_context_log(a->parent->h->context, GRIB_LOG_ERROR,
"grib_accessor_class_bitmap.unpack_bytes : cannot get %s err=%d",self->unusedBits,err);
length-= tlen/8;
memcpy(val,buf + offset,length );
*len = length; *len = length;
return GRIB_ARRAY_TOO_SMALL;
}
if ((err=grib_get_long_internal(a->parent->h, self->unusedBits, &tlen)) != GRIB_SUCCESS) return GRIB_SUCCESS;
grib_context_log(a->parent->h->context, GRIB_LOG_ERROR,
"grib_accessor_class_bitmap.unpack_bytes : cannot get %s err=%d",self->unusedBits,err);
length-= tlen/8;
memcpy(val,buf + offset,length );
*len = length;
return GRIB_SUCCESS;
} }
static void grib_set_bit_on( unsigned char* p, long *bitp){ static void grib_set_bit_on( unsigned char* p, long *bitp){
unsigned char o = 1; unsigned char o = 1;
p += (*bitp >> 3); p += (*bitp >> 3);
o <<= 7-((*bitp)%8); o <<= 7-((*bitp)%8);
*p |= o; *p |= o;
(*bitp)+=1; (*bitp)+=1;
} }

264
src/grib_accessor_class_local_definition.c
View File

@ -9,10 +9,9 @@
*/ */
/************************************** /**************************************
* Enrico Fucile * This pertains to GRIB edition 2
**************************************/ **************************************/
#include "grib_api_internal.h" #include "grib_api_internal.h"
/* /*
This is used by make_class.pl This is used by make_class.pl
@ -150,158 +149,161 @@ static void init_class(grib_accessor_class* c)
static void init(grib_accessor* a,const long l, grib_arguments* c) static void init(grib_accessor* a,const long l, grib_arguments* c)
{ {
grib_accessor_local_definition* self = (grib_accessor_local_definition*)a; grib_accessor_local_definition* self = (grib_accessor_local_definition*)a;
int n = 0; int n = 0;
self->grib2LocalSectionNumber = grib_arguments_get_name(a->parent->h,c,n++);
self->productDefinitionTemplateNumber = grib_arguments_get_name(a->parent->h,c,n++);
self->productDefinitionTemplateNumberInternal = grib_arguments_get_name(a->parent->h,c,n++);
self->type = grib_arguments_get_name(a->parent->h,c,n++);
self->stream = grib_arguments_get_name(a->parent->h,c,n++);
self->class = grib_arguments_get_name(a->parent->h,c,n++);
self->eps = grib_arguments_get_name(a->parent->h,c,n++);
self->stepType = grib_arguments_get_name(a->parent->h,c,n++);
self->derivedForecast = grib_arguments_get_name(a->parent->h,c,n++);
self->grib2LocalSectionNumber = grib_arguments_get_name(a->parent->h,c,n++);
self->productDefinitionTemplateNumber = grib_arguments_get_name(a->parent->h,c,n++);
self->productDefinitionTemplateNumberInternal = grib_arguments_get_name(a->parent->h,c,n++);
self->type = grib_arguments_get_name(a->parent->h,c,n++);
self->stream = grib_arguments_get_name(a->parent->h,c,n++);
self->class = grib_arguments_get_name(a->parent->h,c,n++);
self->eps = grib_arguments_get_name(a->parent->h,c,n++);
self->stepType = grib_arguments_get_name(a->parent->h,c,n++);
self->derivedForecast = grib_arguments_get_name(a->parent->h,c,n++);
} }
static int unpack_long (grib_accessor* a, long* val, size_t *len) static int unpack_long(grib_accessor* a, long* val, size_t *len)
{ {
grib_accessor_local_definition* self = (grib_accessor_local_definition*)a; grib_accessor_local_definition* self = (grib_accessor_local_definition*)a;
return grib_get_long(a->parent->h, self->grib2LocalSectionNumber,val); return grib_get_long(a->parent->h, self->grib2LocalSectionNumber,val);
} }
static int pack_long(grib_accessor* a, const long* val, size_t *len) static int pack_long(grib_accessor* a, const long* val, size_t *len)
{ {
grib_accessor_local_definition* self = (grib_accessor_local_definition*)a; grib_accessor_local_definition* self = (grib_accessor_local_definition*)a;
long productDefinitionTemplateNumber=-1; long productDefinitionTemplateNumber=-1;
long productDefinitionTemplateNumberInternal=-1; long productDefinitionTemplateNumberInternal=-1;
long productDefinitionTemplateNumberNew=-1; long productDefinitionTemplateNumberNew=-1;
long grib2LocalSectionNumber=-1; long grib2LocalSectionNumber=-1;
long type=-1; long type=-1;
long stream=-1; long stream=-1;
long class=-1; long class=-1;
long eps=-1; long eps=-1;
char stepType[15]={0,}; char stepType[15]={0,};
size_t slen=15; size_t slen=15;
int localDefinitionNumber=*val; int localDefinitionNumber=*val;
int isInstant=0; int isInstant=0;
int tooEarly=0; int tooEarly=0;
long derivedForecast=-1; long derivedForecast=-1;
long editionNumber = 0;
if (grib_get_long(a->parent->h, self->productDefinitionTemplateNumber,&productDefinitionTemplateNumber)!=GRIB_SUCCESS) if (grib_get_long(a->parent->h, "editionNumber", &editionNumber)==GRIB_SUCCESS) {
tooEarly=1; Assert(editionNumber == 2);
grib_get_long(a->parent->h, self->productDefinitionTemplateNumberInternal,&productDefinitionTemplateNumberInternal); }
grib_get_long(a->parent->h, self->type,&type);
grib_get_long(a->parent->h, self->stream,&stream);
grib_get_long(a->parent->h, self->class,&class);
grib_get_long(a->parent->h, self->eps,&eps);
grib_get_string(a->parent->h, self->stepType,stepType,&slen);
if (!strcmp(stepType,"instant")) isInstant=1;
grib_get_long(a->parent->h, self->grib2LocalSectionNumber,&grib2LocalSectionNumber);
if (productDefinitionTemplateNumber==1 || productDefinitionTemplateNumber==11) if (grib_get_long(a->parent->h, self->productDefinitionTemplateNumber,&productDefinitionTemplateNumber)!=GRIB_SUCCESS)
eps=1; tooEarly=1;
grib_get_long(a->parent->h, self->productDefinitionTemplateNumberInternal,&productDefinitionTemplateNumberInternal);
grib_get_long(a->parent->h, self->type,&type);
grib_get_long(a->parent->h, self->stream,&stream);
grib_get_long(a->parent->h, self->class,&class);
grib_get_long(a->parent->h, self->eps,&eps);
grib_get_string(a->parent->h, self->stepType,stepType,&slen);
if (!strcmp(stepType,"instant")) isInstant=1;
grib_get_long(a->parent->h, self->grib2LocalSectionNumber,&grib2LocalSectionNumber);
switch (localDefinitionNumber) { if (productDefinitionTemplateNumber==1 || productDefinitionTemplateNumber==11)
case 0: eps=1;
case 300:
productDefinitionTemplateNumberNew=productDefinitionTemplateNumber;
break;
case 500: switch (localDefinitionNumber) {
productDefinitionTemplateNumberNew=0; case 0:
break; case 300:
productDefinitionTemplateNumberNew=productDefinitionTemplateNumber;
break;
case 1: /* MARS labelling */ case 500:
case 36: /* MARS labelling for long window 4Dvar system */ productDefinitionTemplateNumberNew=0;
case 40: /* MARS labeling with domain and model (for LAM) */ break;
if (isInstant) {
/* type=em || type=es */
if (type==17) {
productDefinitionTemplateNumberNew=2;
derivedForecast=0;
} else if (type==18) {
productDefinitionTemplateNumberNew=2;
derivedForecast=4;
/* eps or enda or elda or ewla */
} else if (eps==1 || stream==1030 || stream==1249 || stream==1250) {
productDefinitionTemplateNumberNew=1;
} else {
productDefinitionTemplateNumberNew=0;
}
} else {
/* type=em || type=es */
if (type==17) {
productDefinitionTemplateNumberNew=12;
derivedForecast=0;
} else if (type==18) {
productDefinitionTemplateNumberNew=12;
derivedForecast=4;
/* eps or enda or elda or ewla */
} else if (eps==1 || stream==1030 || stream==1249 || stream==1250) {
productDefinitionTemplateNumberNew=11;
} else {
productDefinitionTemplateNumberNew=8;
}
}
break;
case 15: /* Seasonal forecast data */ case 1: /* MARS labelling */
case 16: /* Seasonal forecast monthly mean data */ case 36: /* MARS labelling for long window 4Dvar system */
case 18: /* Multianalysis ensemble data */ case 40: /* MARS labeling with domain and model (for LAM) */
case 26: /* MARS labelling or ensemble forecast data */ if (isInstant) {
case 30: /* Forecasting Systems with Variable Resolution */ /* type=em || type=es */
if (isInstant) { if (type==17) {
productDefinitionTemplateNumberNew=1; productDefinitionTemplateNumberNew=2;
} else { derivedForecast=0;
productDefinitionTemplateNumberNew=11; } else if (type==18) {
} productDefinitionTemplateNumberNew=2;
break; derivedForecast=4;
/* eps or enda or elda or ewla */
} else if (eps==1 || stream==1030 || stream==1249 || stream==1250) {
productDefinitionTemplateNumberNew=1;
} else {
productDefinitionTemplateNumberNew=0;
}
} else {
/* type=em || type=es */
if (type==17) {
productDefinitionTemplateNumberNew=12;
derivedForecast=0;
} else if (type==18) {
productDefinitionTemplateNumberNew=12;
derivedForecast=4;
/* eps or enda or elda or ewla */
} else if (eps==1 || stream==1030 || stream==1249 || stream==1250) {
productDefinitionTemplateNumberNew=11;
} else {
productDefinitionTemplateNumberNew=8;
}
}
break;
case 7: /* Sensitivity data */ case 15: /* Seasonal forecast data */
case 9: /* Singular vectors and ensemble perturbations */ case 16: /* Seasonal forecast monthly mean data */
case 11: /* Supplementary data used by the analysis */ case 18: /* Multianalysis ensemble data */
case 14: /* Brightness temperature */ case 26: /* MARS labelling or ensemble forecast data */
case 20: /* 4D variational increments */ case 30: /* Forecasting Systems with Variable Resolution */
case 21: /* Sensitive area predictions */ if (isInstant) {
case 23: /* Coupled atmospheric, wave and ocean means */ productDefinitionTemplateNumberNew=1;
case 24: /* Satellite Channel Number Data */ } else {
case 25: productDefinitionTemplateNumberNew=11;
case 28: /* COSMO local area EPS */ }
case 38: /* 4D variational increments for long window 4Dvar system */ break;
case 39: /* 4DVar model errors for long window 4Dvar system */
if (isInstant) {
productDefinitionTemplateNumberNew=0;
} else {
productDefinitionTemplateNumberNew=8;
}
break;
default: case 7: /* Sensitivity data */
grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,"Invalid localDefinitionNumber %d",localDefinitionNumber); case 9: /* Singular vectors and ensemble perturbations */
return GRIB_ENCODING_ERROR; case 11: /* Supplementary data used by the analysis */
break; case 14: /* Brightness temperature */
} case 20: /* 4D variational increments */
case 21: /* Sensitive area predictions */
case 23: /* Coupled atmospheric, wave and ocean means */
case 24: /* Satellite Channel Number Data */
case 25:
case 28: /* COSMO local area EPS */
case 38: /* 4D variational increments for long window 4Dvar system */
case 39: /* 4DVar model errors for long window 4Dvar system */
if (isInstant) {
productDefinitionTemplateNumberNew=0;
} else {
productDefinitionTemplateNumberNew=8;
}
break;
if (productDefinitionTemplateNumber != productDefinitionTemplateNumberNew) { default:
if (tooEarly) grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,"Invalid localDefinitionNumber %d",localDefinitionNumber);
grib_set_long(a->parent->h, self->productDefinitionTemplateNumberInternal,productDefinitionTemplateNumberNew); return GRIB_ENCODING_ERROR;
else break;
grib_set_long(a->parent->h, self->productDefinitionTemplateNumber,productDefinitionTemplateNumberNew); }
}
if (derivedForecast>=0)
grib_set_long(a->parent->h, self->derivedForecast,derivedForecast);
grib_set_long(a->parent->h, self->grib2LocalSectionNumber,*val); if (productDefinitionTemplateNumber != productDefinitionTemplateNumberNew) {
if (tooEarly)
grib_set_long(a->parent->h, self->productDefinitionTemplateNumberInternal,productDefinitionTemplateNumberNew);
else
grib_set_long(a->parent->h, self->productDefinitionTemplateNumber,productDefinitionTemplateNumberNew);
}
if (derivedForecast>=0)
grib_set_long(a->parent->h, self->derivedForecast,derivedForecast);
return 0; grib_set_long(a->parent->h, self->grib2LocalSectionNumber,*val);
return 0;
} }
static long value_count(grib_accessor* a) static long value_count(grib_accessor* a)
{ {
return 1; return 1;
} }

314
src/grib_dumper_class_json.c
View File

@ -84,7 +84,7 @@ grib_dumper_class* grib_dumper_class_json = &_grib_dumper_class_json;
/* END_CLASS_IMP */ /* END_CLASS_IMP */
GRIB_INLINE static int grib_inline_strcmp(const char* a,const char* b) { GRIB_INLINE static int grib_inline_strcmp(const char* a,const char* b) {
if (*a != *b) return 1; if (*a != *b) return 1;
while((*a!=0 && *b!=0) && *(a) == *(b) ) {a++;b++;} while((*a!=0 && *b!=0) && *(a) == *(b) ) {a++;b++;}
return (*a==0 && *b==0) ? 0 : 1; return (*a==0 && *b==0) ? 0 : 1;
} }
@ -95,147 +95,147 @@ static void init_class (grib_dumper_class* c){}
static int init(grib_dumper* d) static int init(grib_dumper* d)
{ {
grib_dumper_json *self = (grib_dumper_json*)d; grib_dumper_json *self = (grib_dumper_json*)d;
self->section_offset=0; self->section_offset=0;
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static int destroy (grib_dumper* d) static int destroy (grib_dumper* d)
{ {
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static void dump_values(grib_dumper* d,grib_accessor* a) static void dump_values(grib_dumper* d,grib_accessor* a)
{ {
grib_dumper_json *self = (grib_dumper_json*)d; grib_dumper_json *self = (grib_dumper_json*)d;
double value; size_t size = 1; double value; size_t size = 1;
double *values=NULL; double *values=NULL;
int err = 0; int err = 0;
int i,tab; int i,tab;
long more=0; long more=0;
long count=0; long count=0;
int mydepth=depth+2; int mydepth=depth+2;
double missing_value = 9999; double missing_value = 9999;
count=grib_value_count(a); count=grib_value_count(a);
size=count; size=count;
if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0) if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0)
return; return;
if (size>1) { if (size>1) {
values=grib_context_malloc_clear(a->parent->h->context,sizeof(double)*size); values=grib_context_malloc_clear(a->parent->h->context,sizeof(double)*size);
err=grib_unpack_double(a,values,&size); err=grib_unpack_double(a,values,&size);
} else { } else {
err=grib_unpack_double(a,&value,&size); err=grib_unpack_double(a,&value,&size);
} }
if(!(d->option_flags & GRIB_DUMP_FLAG_ALL_DATA) && size > 3) { if(!(d->option_flags & GRIB_DUMP_FLAG_ALL_DATA) && size > 3) {
more = size - 3; more = size - 3;
size = 3; size = 3;
} else more=0; } else more=0;
if (!self->begin) fprintf(self->dumper.out,",\n"); if (!self->begin) fprintf(self->dumper.out,",\n");
else self->begin=0; else self->begin=0;
if (size>1) { if (size>1) {
int cols=4; int cols=4;
int count=0; int count=0;
int lens=strlen(a->name); int lens=strlen(a->name);
fprintf(self->dumper.out,"%-*s",mydepth," "); fprintf(self->dumper.out,"%-*s",mydepth," ");
fprintf(self->dumper.out,"\"%s\" : [ ",a->name); fprintf(self->dumper.out,"\"%s\" : [ ",a->name);
tab=lens+mydepth+7; tab=lens+mydepth+7;
err = grib_get_double(a->parent->h, "missingValue", &missing_value); err = grib_get_double(a->parent->h, "missingValue", &missing_value);
for (i=0; i<size-1; ++i) { for (i=0; i<size-1; ++i) {
if (count>cols || i==0) {fprintf(self->dumper.out,"\n%-*s",tab," ");count=0;} if (count>cols || i==0) {fprintf(self->dumper.out,"\n%-*s",tab," ");count=0;}
if (values[i] == missing_value) if (values[i] == missing_value)
fprintf(self->dumper.out, "%s, ", "null"); fprintf(self->dumper.out, "%s, ", "null");
else else
fprintf(self->dumper.out,"%g, ",values[i]); fprintf(self->dumper.out, "%g, ", values[i]);
count++; count++;
} }
if (count>cols) fprintf(self->dumper.out,"\n%-*s",tab," "); if (count>cols) fprintf(self->dumper.out,"\n%-*s",tab," ");
if (values[i] == missing_value) if (values[i] == missing_value)
fprintf(self->dumper.out, "%s ","null"); fprintf(self->dumper.out, "%s ","null");
else else
fprintf(self->dumper.out,"%g ",values[i]); fprintf(self->dumper.out, "%g ",values[i]);
if (more) if (more)
fprintf(self->dumper.out,"\n%-*s... %ld more values",tab," ",more); fprintf(self->dumper.out, "\n%-*s... %ld more values", tab, " ", more);
tab=lens+mydepth+5; tab=lens+mydepth+5;
fprintf(self->dumper.out,"\n%-*s] ",tab," "); fprintf(self->dumper.out,"\n%-*s] ",tab," ");
grib_context_free(a->parent->h->context,values); grib_context_free(a->parent->h->context,values);
} else { } else {
fprintf(self->dumper.out,"%-*s",mydepth," "); fprintf(self->dumper.out,"%-*s",mydepth," ");
if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) ) if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) )
fprintf(self->dumper.out,"\"%s\" : null",a->name); fprintf(self->dumper.out,"\"%s\" : null",a->name);
else else
fprintf(self->dumper.out,"\"%s\" : %g",a->name,value); fprintf(self->dumper.out,"\"%s\" : %g",a->name,value);
} }
} }
static void dump_long(grib_dumper* d,grib_accessor* a,const char* comment) static void dump_long(grib_dumper* d,grib_accessor* a,const char* comment)
{ {
grib_dumper_json *self = (grib_dumper_json*)d; grib_dumper_json *self = (grib_dumper_json*)d;
long value; size_t size = 1; long value; size_t size = 1;
long *values=NULL; long *values=NULL;
int err = 0; int err = 0;
int i,tab; int i,tab;
long count=0; long count=0;
long more=0; long more=0;
int mydepth=depth+2; int mydepth=depth+2;
count = grib_value_count(a); count = grib_value_count(a);
size=count; size=count;
if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0) if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0)
return; return;
if (size>1) { if (size>1) {
values=grib_context_malloc_clear(a->parent->h->context,sizeof(long)*size); values=grib_context_malloc_clear(a->parent->h->context,sizeof(long)*size);
err=grib_unpack_long(a,values,&size); err=grib_unpack_long(a,values,&size);
} else { } else {
err=grib_unpack_long(a,&value,&size); err=grib_unpack_long(a,&value,&size);
}
if(!(d->option_flags & GRIB_DUMP_FLAG_ALL_DATA) && size > 3) {
more = size - 3;
size = 3;
} else more=0;
if (!self->begin) fprintf(self->dumper.out,",\n");
else self->begin=0;
if (size>1) {
int cols=9;
int count=0;
int lens=strlen(a->name);
fprintf(self->dumper.out,"%-*s",mydepth," ");
fprintf(self->dumper.out,"\"%s\" : [ ",a->name);
tab=lens+mydepth+7;
for (i=0;i<size-1;i++) {
if (count>cols || i==0) {fprintf(self->dumper.out,"\n%-*s",tab," ");count=0;}
fprintf(self->dumper.out,"%ld, ",values[i]);
count++;
} }
if (count>cols) fprintf(self->dumper.out,"\n%-*s",tab," ");
fprintf(self->dumper.out,"%ld ",values[i]);
if (more)
fprintf(self->dumper.out,"\n%-*s... %ld more values",tab," ",more);
tab=lens+mydepth+5; if(!(d->option_flags & GRIB_DUMP_FLAG_ALL_DATA) && size > 3) {
fprintf(self->dumper.out,"\n%-*s] ",tab," "); more = size - 3;
grib_context_free(a->parent->h->context,values); size = 3;
} else { } else more=0;
fprintf(self->dumper.out,"%-*s",mydepth," ");
if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) ) if (!self->begin) fprintf(self->dumper.out,",\n");
fprintf(self->dumper.out,"\"%s\" : null",a->name); else self->begin=0;
else
fprintf(self->dumper.out,"\"%s\" : %ld",a->name,value); if (size>1) {
} int cols=9;
int count=0;
int lens=strlen(a->name);
fprintf(self->dumper.out,"%-*s",mydepth," ");
fprintf(self->dumper.out,"\"%s\" : [ ",a->name);
tab=lens+mydepth+7;
for (i=0;i<size-1;i++) {
if (count>cols || i==0) {fprintf(self->dumper.out,"\n%-*s",tab," ");count=0;}
fprintf(self->dumper.out,"%ld, ",values[i]);
count++;
}
if (count>cols) fprintf(self->dumper.out,"\n%-*s",tab," ");
fprintf(self->dumper.out,"%ld ",values[i]);
if (more)
fprintf(self->dumper.out,"\n%-*s... %ld more values",tab," ",more);
tab=lens+mydepth+5;
fprintf(self->dumper.out,"\n%-*s] ",tab," ");
grib_context_free(a->parent->h->context,values);
} else {
fprintf(self->dumper.out,"%-*s",mydepth," ");
if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) )
fprintf(self->dumper.out,"\"%s\" : null",a->name);
else
fprintf(self->dumper.out,"\"%s\" : %ld",a->name,value);
}
} }
@ -245,62 +245,62 @@ static void dump_bits(grib_dumper* d,grib_accessor* a,const char* comment)
static void dump_double(grib_dumper* d,grib_accessor* a,const char* comment) static void dump_double(grib_dumper* d,grib_accessor* a,const char* comment)
{ {
grib_dumper_json *self = (grib_dumper_json*)d; grib_dumper_json *self = (grib_dumper_json*)d;
double value; size_t size = 1; double value; size_t size = 1;
int mydepth=depth+2; int mydepth=depth+2;
grib_unpack_double(a,&value,&size); grib_unpack_double(a,&value,&size);
if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0) if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0)
return; return;
if (!self->begin) fprintf(self->dumper.out,",\n"); if (!self->begin) fprintf(self->dumper.out,",\n");
else self->begin=0; else self->begin=0;
fprintf(self->dumper.out,"%-*s",mydepth," "); fprintf(self->dumper.out,"%-*s",mydepth," ");
if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) ) if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) )
fprintf(self->dumper.out,"\"%s\" : null",a->name); fprintf(self->dumper.out,"\"%s\" : null",a->name);
else else
fprintf(self->dumper.out,"\"%s\" : %g",a->name,value); fprintf(self->dumper.out,"\"%s\" : %g",a->name,value);
} }
static void dump_string(grib_dumper* d,grib_accessor* a,const char* comment) static void dump_string(grib_dumper* d,grib_accessor* a,const char* comment)
{ {
grib_dumper_json *self = (grib_dumper_json*)d; grib_dumper_json *self = (grib_dumper_json*)d;
char *value=NULL; char *value=NULL;
char *p = NULL; char *p = NULL;
size_t size = 0; size_t size = 0;
grib_context* c=NULL; grib_context* c=NULL;
int err = grib_get_string_length(a->parent->h,a->name,&size); int err = grib_get_string_length(a->parent->h,a->name,&size);
int mydepth=depth+2; int mydepth=depth+2;
c=a->parent->h->context; c=a->parent->h->context;
if (size==0) return; if (size==0) return;
if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0) if ( (a->flags & GRIB_ACCESSOR_FLAG_DUMP) == 0)
return; return;
value=grib_context_malloc_clear(c,size); value=grib_context_malloc_clear(c,size);
if (!value) { if (!value) {
grib_context_log(c,GRIB_LOG_FATAL,"unable to allocate %d bytes",(int)size); grib_context_log(c,GRIB_LOG_FATAL,"unable to allocate %d bytes",(int)size);
return; return;
} }
if (!self->begin) fprintf(self->dumper.out,",\n"); if (!self->begin) fprintf(self->dumper.out,",\n");
else self->begin=0; else self->begin=0;
err = grib_unpack_string(a,value,&size); err = grib_unpack_string(a,value,&size);
p=value; p=value;
while(*p) { if(!isprint(*p)) *p = '.'; p++; } while(*p) { if(!isprint(*p)) *p = '.'; p++; }
fprintf(self->dumper.out,"%-*s",mydepth," "); fprintf(self->dumper.out,"%-*s",mydepth," ");
if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) ) if( ((a->flags & GRIB_ACCESSOR_FLAG_CAN_BE_MISSING) != 0) && grib_is_missing_internal(a) )
fprintf(self->dumper.out,"\"%s\" : null",a->name); fprintf(self->dumper.out,"\"%s\" : null",a->name);
else else
fprintf(self->dumper.out,"\"%s\" : \"%s\"",a->name,value); fprintf(self->dumper.out,"\"%s\" : \"%s\"",a->name,value);
grib_context_free(c,value); grib_context_free(c,value);
} }
static void dump_bytes(grib_dumper* d,grib_accessor* a,const char* comment) static void dump_bytes(grib_dumper* d,grib_accessor* a,const char* comment)
@ -313,15 +313,15 @@ static void dump_label(grib_dumper* d,grib_accessor* a,const char* comment)
static void dump_section(grib_dumper* d,grib_accessor* a,grib_block_of_accessors* block) static void dump_section(grib_dumper* d,grib_accessor* a,grib_block_of_accessors* block)
{ {
grib_dumper_json *self = (grib_dumper_json*)d; grib_dumper_json *self = (grib_dumper_json*)d;
if ( !grib_inline_strcmp(a->name,"GRIB") ) { if ( !grib_inline_strcmp(a->name,"GRIB") ) {
fprintf(self->dumper.out,"{\n"); fprintf(self->dumper.out,"{\n");
self->begin=1; self->begin=1;
grib_dump_accessors_block(d,block); grib_dump_accessors_block(d,block);
fprintf(self->dumper.out,"\n}\n"); fprintf(self->dumper.out,"\n}\n");
} }
else { else {
grib_dump_accessors_block(d,block); grib_dump_accessors_block(d,block);
} }
} }

2
src/grib_errors.c.in
View File

@ -35,7 +35,7 @@ void grib_check(const char* call,const char* file,int line,int e,const char* ms
} }
} }
void grib_fail(const char* expr,const char* file,int line) { void grib_fail(const char* expr,const char* file,int line) {
fprintf(stderr,"%s at line %d: assertion failure Assert(%s)\n",file,line,expr); fprintf(stderr,"%s at line %d: assertion failure Assert(%s)\n",file,line,expr);
abort(); abort();
} }

280
src/grib_ibmfloat.c
View File

@ -15,212 +15,212 @@ static void init_ibm_table();
typedef struct ibm_table_t ibm_table_t; typedef struct ibm_table_t ibm_table_t;
struct ibm_table_t { struct ibm_table_t {
int inited; int inited;
double e[128]; double e[128];
double v[128]; double v[128];
double vmin; double vmin;
double vmax; double vmax;
}; };
static ibm_table_t ibm_table={ 0, {0,}, {0,}, 0, 0 }; static ibm_table_t ibm_table={ 0, {0,}, {0,}, 0, 0 };
static void init_ibm_table() { static void init_ibm_table() {
if (!ibm_table.inited) { if (!ibm_table.inited) {
unsigned long i; unsigned long i;
unsigned long mmin=0x100000; unsigned long mmin=0x100000;
unsigned long mmax=0xffffff; unsigned long mmax=0xffffff;
double e=1; double e=1;
for (i=1; i<=57;i++) { e*=16; ibm_table.e[i+70]=e;ibm_table.v[i+70]=e*mmin;} for (i=1; i<=57;i++) { e*=16; ibm_table.e[i+70]=e;ibm_table.v[i+70]=e*mmin;}
ibm_table.e[70]=1; ibm_table.e[70]=1;
ibm_table.v[70]=mmin; ibm_table.v[70]=mmin;
e=1; e=1;
for (i=1; i<=70;i++) { e/=16; ibm_table.e[-i+70]=e;ibm_table.v[-i+70]=e*mmin;} for (i=1; i<=70;i++) { e/=16; ibm_table.e[-i+70]=e;ibm_table.v[-i+70]=e*mmin;}
ibm_table.vmin=ibm_table.v[0]; ibm_table.vmin=ibm_table.v[0];
ibm_table.vmax=ibm_table.e[127]*mmax; ibm_table.vmax=ibm_table.e[127]*mmax;
ibm_table.inited=1; ibm_table.inited=1;
/*for (i=0;i<128;i++) printf("++++ ibm_table.v[%d]=%g\n",i,ibm_table.v[i]);*/ /*for (i=0;i<128;i++) printf("++++ ibm_table.v[%d]=%g\n",i,ibm_table.v[i]);*/
} }
} }
static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j) static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j)
{ {
/*These routine works only on ascending ordered arrays*/ /*These routine works only on ascending ordered arrays*/
unsigned long ju,jm,jl; unsigned long ju,jm,jl;
jl=0; jl=0;
ju=n; ju=n;
while (ju-jl > 1) { while (ju-jl > 1) {
jm=(ju+jl) >> 1; jm=(ju+jl) >> 1;
/* printf("jl=%lu jm=%lu ju=%lu\n",jl,jm,ju); */ /* printf("jl=%lu jm=%lu ju=%lu\n",jl,jm,ju); */
/* printf("xx[jl]=%.10e xx[jm]=%.10e xx[ju]=%.10e\n",xx[jl],xx[jm],xx[ju]); */ /* printf("xx[jl]=%.10e xx[jm]=%.10e xx[ju]=%.10e\n",xx[jl],xx[jm],xx[ju]); */
if (x >= xx[jm]) jl=jm; if (x >= xx[jm]) jl=jm;
else ju=jm; else ju=jm;
} }
*j=jl; *j=jl;
} }
unsigned long grib_ibm_to_long(double x) { unsigned long grib_ibm_to_long(double x) {
unsigned long s = 0; unsigned long s = 0;
unsigned long mmax = 0xffffff; unsigned long mmax = 0xffffff;
unsigned long mmin = 0x800000; unsigned long mmin = 0x800000;
unsigned long m = mmax; unsigned long m = mmax;
unsigned long e=0; unsigned long e=0;
double rmmax=mmax+0.5; double rmmax=mmax+0.5;
if (!ibm_table.inited) init_ibm_table(); if (!ibm_table.inited) init_ibm_table();
/* printf("\ngrib_ibm_to_long: x=%.20e\n",x); */ /* printf("\ngrib_ibm_to_long: x=%.20e\n",x); */
if (x < 0) { s = 1; x = -x; } if (x < 0) { s = 1; x = -x; }
/* Underflow */ /* Underflow */
if (x < ibm_table.vmin) { if (x < ibm_table.vmin) {
/*printf("grib_ibm_to_long: (x < ibm_table.vmin) x=%.20e vmin=%.20e v=0x%lX\n",x,ibm_table.vmin,(s<<31));*/ /*printf("grib_ibm_to_long: (x < ibm_table.vmin) x=%.20e vmin=%.20e v=0x%lX\n",x,ibm_table.vmin,(s<<31));*/
return (s << 31); return (s << 31);
} }
/* Overflow */ /* Overflow */
if (x > ibm_table.vmax) { if (x > ibm_table.vmax) {
fprintf(stderr, "grib_ibm_to_long: Number is too large: x=%.20e > xmax=%.20e\n", x, ibm_table.vmax); fprintf(stderr, "grib_ibm_to_long: Number is too large: x=%.20e > xmax=%.20e\n", x, ibm_table.vmax);
Assert(0); Assert(0);
return 0; return 0;
} }
binary_search(ibm_table.v, 127, x, &e); binary_search(ibm_table.v, 127, x, &e);
/* printf("grib_ibm_to_long: e=%ld\n",e); */ /* printf("grib_ibm_to_long: e=%ld\n",e); */
x/=ibm_table.e[e]; x/=ibm_table.e[e];
/* printf("grib_ibm_to_long: x=%.20e\n",x); */ /* printf("grib_ibm_to_long: x=%.20e\n",x); */
while(x < mmin ) { x *= 16; e--; while(x < mmin ) { x *= 16; e--;
/* printf("grib_ibm_to_long (e--): x=%.20e e=%ld \n",x,e); */ /* printf("grib_ibm_to_long (e--): x=%.20e e=%ld \n",x,e); */
} }
while(x > rmmax ) { x /= 16; e++; while(x > rmmax ) { x /= 16; e++;
/* printf("grib_ibm_to_long (e++): x=%.20e e=%ld \n",x,e); */ /* printf("grib_ibm_to_long (e++): x=%.20e e=%ld \n",x,e); */
} }
m=x+0.5; m=x+0.5;
/* printf("grib_ibm_to_long: m=0x%lX (%lu) x=%.10e \n",m,m,x ); */ /* printf("grib_ibm_to_long: m=0x%lX (%lu) x=%.10e \n",m,m,x ); */
if ( m > mmax ) { e++; m=0x800000; if ( m > mmax ) { e++; m=0x800000;
/* printf("grib_ibm_to_long: ( m > mmax ) m=0x%lX (%lu) x=%.10e \n",m,m,x ); */ /* printf("grib_ibm_to_long: ( m > mmax ) m=0x%lX (%lu) x=%.10e \n",m,m,x ); */
} }
/* printf("grib_ibm_to_long: s=%lu c=%lu (0x%lX) m=%lu (0x%lX)\n",s,e,e,m,m ); */ /* printf("grib_ibm_to_long: s=%lu c=%lu (0x%lX) m=%lu (0x%lX)\n",s,e,e,m,m ); */
return (s << 31) | ( e << 24 ) | m; return (s << 31) | ( e << 24 ) | m;
} }
double grib_ibmfloat_error(double x) { double grib_ibmfloat_error(double x) {
unsigned long e=0; unsigned long e=0;
if (!ibm_table.inited) init_ibm_table();
if (x < 0) x = -x; if (!ibm_table.inited) init_ibm_table();
/* Underflow */ if (x < 0) x = -x;
if (x <= ibm_table.vmin) return ibm_table.vmin;
/* Overflow */ /* Underflow */
if (x > ibm_table.vmax) { if (x <= ibm_table.vmin) return ibm_table.vmin;
fprintf(stderr, "grib_ibmfloat_error: Number is too large: x=%.20e > xmax=%.20e\n", x, ibm_table.vmax);
Assert(0);
return 0;
}
binary_search(ibm_table.v, 127, x, &e); /* Overflow */
if (x > ibm_table.vmax) {
return ibm_table.e[e] ; fprintf(stderr, "grib_ibmfloat_error: Number is too large: x=%.20e > xmax=%.20e\n", x, ibm_table.vmax);
Assert(0);
return 0;
}
binary_search(ibm_table.v, 127, x, &e);
return ibm_table.e[e] ;
} }
double grib_long_to_ibm(unsigned long x){ double grib_long_to_ibm(unsigned long x){
unsigned long s = x & 0x80000000; unsigned long s = x & 0x80000000;
unsigned long c = (x & 0x7f000000) >> 24; unsigned long c = (x & 0x7f000000) >> 24;
unsigned long m = (x & 0x00ffffff); unsigned long m = (x & 0x00ffffff);
double val = m; double val = m;
if (!ibm_table.inited) init_ibm_table(); if (!ibm_table.inited) init_ibm_table();
/*if(x == 0) return 0;*/ /*if(x == 0) return 0;*/
if (c==0 && m <= 1 ) return 0; if (c==0 && m <= 1 ) return 0;
val*=ibm_table.e[c]; val*=ibm_table.e[c];
if(s) val = -val; if(s) val = -val;
return val; return val;
} }
double grib_ibm_table_e(unsigned long e) { double grib_ibm_table_e(unsigned long e) {
if (!ibm_table.inited) init_ibm_table(); if (!ibm_table.inited) init_ibm_table();
return ibm_table.e[e]; return ibm_table.e[e];
} }
double grib_ibm_table_v(unsigned long e) { double grib_ibm_table_v(unsigned long e) {
if (!ibm_table.inited) init_ibm_table(); if (!ibm_table.inited) init_ibm_table();
return ibm_table.v[e]; return ibm_table.v[e];
} }
unsigned long grib_ibm_nearest_smaller_to_long(double x) unsigned long grib_ibm_nearest_smaller_to_long(double x)
{ {
unsigned long l; unsigned long l;
unsigned long e; unsigned long e;
unsigned long m ; unsigned long m ;
unsigned long s; unsigned long s;
unsigned long mmin = 0x100000; unsigned long mmin = 0x100000;
double y, eps=0; double y, eps=0;
if(x == 0) return 0; if(x == 0) return 0;
if (!ibm_table.inited) init_ibm_table(); if (!ibm_table.inited) init_ibm_table();
l=grib_ibm_to_long(x); l=grib_ibm_to_long(x);
y=grib_long_to_ibm(l); y=grib_long_to_ibm(l);
if ( x < y ) { if ( x < y ) {
if ( x < 0 && -x < ibm_table.vmin ) { if ( x < 0 && -x < ibm_table.vmin ) {
l=0x80100000; l=0x80100000;
} else { } else {
e = (l & 0x7f000000) >> 24; e = (l & 0x7f000000) >> 24;
m = (l & 0x00ffffff); m = (l & 0x00ffffff);
s = l & 0x80000000; s = l & 0x80000000;
if ( m == mmin ) { if ( m == mmin ) {
/* printf("grib_ibm_nearest_smaller_to_long: m == mmin (0x%lX) e=%lu\n",m,e); */ /* printf("grib_ibm_nearest_smaller_to_long: m == mmin (0x%lX) e=%lu\n",m,e); */
e = s ? e : e-1; e = s ? e : e-1;
/*if (e<0) e=0; this condition is always FALSE -- 'e' is unsigned long */ /*if (e<0) e=0; this condition is always FALSE -- 'e' is unsigned long */
if (e>127) e=127; if (e>127) e=127;
/* printf("grib_ibm_nearest_smaller_to_long: e=%lu \n",e); */ /* printf("grib_ibm_nearest_smaller_to_long: e=%lu \n",e); */
} }
eps=ibm_table.e[e]; eps=ibm_table.e[e];
/* printf("grib_ibm_nearest_smaller_to_long: x<y\n"); */ /* printf("grib_ibm_nearest_smaller_to_long: x<y\n"); */
l=grib_ibm_to_long(y-eps); l=grib_ibm_to_long(y-eps);
/* printf("grib_ibm_nearest_smaller_to_long: grib_ibm_to_long(y-eps)=0x%lX y=%.10e eps=%.10e x=%.10e\n",l,y,eps,x); */ /* printf("grib_ibm_nearest_smaller_to_long: grib_ibm_to_long(y-eps)=0x%lX y=%.10e eps=%.10e x=%.10e\n",l,y,eps,x); */
} }
} }
if (x<grib_long_to_ibm(l)) { if (x<grib_long_to_ibm(l)) {
l=grib_ibm_to_long(x-eps); l=grib_ibm_to_long(x-eps);
if (x<grib_long_to_ibm(l)) { if (x<grib_long_to_ibm(l)) {
printf("grib_ibm_nearest_smaller_to_long: x=%.20e grib_long_to_ibm(0x%lX)=%.20e\n",x,l,grib_long_to_ibm(l)); printf("grib_ibm_nearest_smaller_to_long: x=%.20e grib_long_to_ibm(0x%lX)=%.20e\n",x,l,grib_long_to_ibm(l));
Assert(x >= grib_long_to_ibm(l)); Assert(x >= grib_long_to_ibm(l));
} }
} }
return l; return l;
} }
int grib_nearest_smaller_ibm_float(double a,double* ret) { int grib_nearest_smaller_ibm_float(double a,double* ret) {
unsigned long l=0; unsigned long l=0;
if (!ibm_table.inited) init_ibm_table(); if (!ibm_table.inited) init_ibm_table();
if (a>ibm_table.vmax) return GRIB_INTERNAL_ERROR; if (a>ibm_table.vmax) return GRIB_INTERNAL_ERROR;
l=grib_ibm_nearest_smaller_to_long(a); l=grib_ibm_nearest_smaller_to_long(a);
*ret=grib_long_to_ibm(l); *ret=grib_long_to_ibm(l);
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }

594
src/grib_ieeefloat.c
View File

@ -21,221 +21,221 @@ static void init_ieee_table();
typedef struct ieee_table_t ieee_table_t; typedef struct ieee_table_t ieee_table_t;
struct ieee_table_t { struct ieee_table_t {
int inited; int inited;
double e[255]; double e[255];
double v[255]; double v[255];
double vmin; double vmin;
double vmax; double vmax;
}; };
static ieee_table_t ieee_table={ 0,{0,},{0,}, 0, 0 }; static ieee_table_t ieee_table={ 0,{0,},{0,}, 0, 0 };
static void init_ieee_table() { static void init_ieee_table() {
if (!ieee_table.inited) { if (!ieee_table.inited) {
unsigned long i; unsigned long i;
unsigned long mmin=0x800000; unsigned long mmin=0x800000;
unsigned long mmax=0xffffff; unsigned long mmax=0xffffff;
double e=1; double e=1;
for (i=1; i<=104;i++) { for (i=1; i<=104;i++) {
e*=2; e*=2;
ieee_table.e[i+150]=e; ieee_table.e[i+150]=e;
ieee_table.v[i+150]=e*mmin; ieee_table.v[i+150]=e*mmin;
}
ieee_table.e[150]=1;
ieee_table.v[150]=mmin;
e=1;
for (i=1; i<150;i++) {
e/=2;
ieee_table.e[-i+150]=e;
ieee_table.v[-i+150]=e*mmin;
}
ieee_table.vmin=ieee_table.v[1];
ieee_table.vmax=ieee_table.e[254]*mmax;
ieee_table.inited=1;
/*for (i=0;i<128;i++) printf("++++ ieee_table.v[%d]=%g\n",i,ieee_table.v[i]);*/
} }
ieee_table.e[150]=1;
ieee_table.v[150]=mmin;
e=1;
for (i=1; i<150;i++) {
e/=2;
ieee_table.e[-i+150]=e;
ieee_table.v[-i+150]=e*mmin;
}
ieee_table.vmin=ieee_table.v[1];
ieee_table.vmax=ieee_table.e[254]*mmax;
ieee_table.inited=1;
/*for (i=0;i<128;i++) printf("++++ ieee_table.v[%d]=%g\n",i,ieee_table.v[i]);*/
}
} }
static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j) static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j)
{ {
/*These routine works only on ascending ordered arrays*/ /*These routine works only on ascending ordered arrays*/
unsigned long ju,jm,jl; unsigned long ju,jm,jl;
jl=0; jl=0;
ju=n; ju=n;
while (ju-jl > 1) { while (ju-jl > 1) {
jm=(ju+jl) >> 1; jm=(ju+jl) >> 1;
/* printf("jl=%lu jm=%lu ju=%lu\n",jl,jm,ju); */ /* printf("jl=%lu jm=%lu ju=%lu\n",jl,jm,ju); */
/* printf("xx[jl]=%.10e xx[jm]=%.10e xx[ju]=%.10e\n",xx[jl],xx[jm],xx[ju]); */ /* printf("xx[jl]=%.10e xx[jm]=%.10e xx[ju]=%.10e\n",xx[jl],xx[jm],xx[ju]); */
if (x >= xx[jm]) jl=jm; if (x >= xx[jm]) jl=jm;
else ju=jm; else ju=jm;
} }
*j=jl; *j=jl;
} }
double grib_ieee_table_e(unsigned long e) { double grib_ieee_table_e(unsigned long e) {
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
return ieee_table.e[e]; return ieee_table.e[e];
} }
double grib_ieee_table_v(unsigned long e) { double grib_ieee_table_v(unsigned long e) {
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
return ieee_table.v[e]; return ieee_table.v[e];
} }
unsigned long grib_ieee_to_long(double x) { unsigned long grib_ieee_to_long(double x) {
unsigned long s = 0; unsigned long s = 0;
unsigned long mmax = 0xffffff; unsigned long mmax = 0xffffff;
unsigned long mmin = 0x800000; unsigned long mmin = 0x800000;
unsigned long m = mmax; unsigned long m = mmax;
unsigned long e=0; unsigned long e=0;
double rmmax=mmax+0.5; double rmmax=mmax+0.5;
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
/* printf("\ngrib_ieee_to_long: x=%.20e\n",x); */ /* printf("\ngrib_ieee_to_long: x=%.20e\n",x); */
if (x < 0) { s = 1; x = -x; } if (x < 0) { s = 1; x = -x; }
/* Underflow */ /* Underflow */
if (x < ieee_table.vmin) { if (x < ieee_table.vmin) {
/*printf("grib_ieee_to_long: (x < ieee_table.vmin) x=%.20e vmin=%.20e v=0x%lX\n",x,ieee_table.vmin,(s<<31));*/ /*printf("grib_ieee_to_long: (x < ieee_table.vmin) x=%.20e vmin=%.20e v=0x%lX\n",x,ieee_table.vmin,(s<<31));*/
return (s << 31); return (s << 31);
} }
/* Overflow */ /* Overflow */
if (x > ieee_table.vmax) { if (x > ieee_table.vmax) {
fprintf(stderr, "grib_ieee_to_long: Number is too large: x=%.20e > xmax=%.20e\n", x, ieee_table.vmax); fprintf(stderr, "grib_ieee_to_long: Number is too large: x=%.20e > xmax=%.20e\n", x, ieee_table.vmax);
Assert(0); Assert(0);
return 0; return 0;
} }
binary_search(ieee_table.v, 254, x, &e); binary_search(ieee_table.v, 254, x, &e);
/* printf("grib_ieee_to_long: e=%ld\n",e); */ /* printf("grib_ieee_to_long: e=%ld\n",e); */
x/=ieee_table.e[e]; x/=ieee_table.e[e];
/* printf("grib_ieee_to_long: x=%.20e\n",x); */ /* printf("grib_ieee_to_long: x=%.20e\n",x); */
while(x < mmin ) { x *= 2; e--; while(x < mmin ) { x *= 2; e--;
/* printf("grib_ieee_to_long (e--): x=%.20e e=%ld \n",x,e); */ /* printf("grib_ieee_to_long (e--): x=%.20e e=%ld \n",x,e); */
} }
while(x > rmmax ) { x /= 2; e++; while(x > rmmax ) { x /= 2; e++;
/* printf("grib_ieee_to_long (e++): x=%.20e e=%ld \n",x,e); */ /* printf("grib_ieee_to_long (e++): x=%.20e e=%ld \n",x,e); */
} }
m=x+0.5; m=x+0.5;
/* printf("grib_ieee_to_long: m=0x%lX (%lu) x=%.10e \n",m,m,x ); */ /* printf("grib_ieee_to_long: m=0x%lX (%lu) x=%.10e \n",m,m,x ); */
if ( m > mmax ) { e++; m=0x800000; if ( m > mmax ) { e++; m=0x800000;
/* printf("grib_ieee_to_long: ( m > mmax ) m=0x%lX (%lu) x=%.10e \n",m,m,x ); */ /* printf("grib_ieee_to_long: ( m > mmax ) m=0x%lX (%lu) x=%.10e \n",m,m,x ); */
} }
/* printf("grib_ieee_to_long: s=%lu c=%lu (0x%lX) m=%lu (0x%lX)\n",s,e,e,m,m ); */ /* printf("grib_ieee_to_long: s=%lu c=%lu (0x%lX) m=%lu (0x%lX)\n",s,e,e,m,m ); */
return (s << 31) | ( e << 23 ) | ( m & 0x7fffff ); return (s << 31) | ( e << 23 ) | ( m & 0x7fffff );
} }
double grib_ieeefloat_error(double x) { double grib_ieeefloat_error(double x) {
unsigned long e=0; unsigned long e=0;
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
if (x < 0) x = -x; if (x < 0) x = -x;
/* Underflow */ /* Underflow */
if (x < ieee_table.vmin) return ieee_table.vmin; if (x < ieee_table.vmin) return ieee_table.vmin;
/* Overflow */ /* Overflow */
if (x > ieee_table.vmax) { if (x > ieee_table.vmax) {
fprintf(stderr, "grib_ieeefloat_error: Number is too large: x=%.20e > xmax=%.20e\n", x, ieee_table.vmax); fprintf(stderr, "grib_ieeefloat_error: Number is too large: x=%.20e > xmax=%.20e\n", x, ieee_table.vmax);
Assert(0); Assert(0);
return 0; return 0;
} }
binary_search(ieee_table.v, 254, x, &e); binary_search(ieee_table.v, 254, x, &e);
return ieee_table.e[e] ; return ieee_table.e[e] ;
} }
double grib_long_to_ieee(unsigned long x){ double grib_long_to_ieee(unsigned long x){
unsigned long s = x & 0x80000000; unsigned long s = x & 0x80000000;
unsigned long c = (x & 0x7f800000) >> 23; unsigned long c = (x & 0x7f800000) >> 23;
unsigned long m = (x & 0x007fffff); unsigned long m = (x & 0x007fffff);
double val; double val;
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
if (c == 0 && m==0) return 0; if (c == 0 && m==0) return 0;
if (c == 0) { if (c == 0) {
m |= 0x800000; m |= 0x800000;
c=1; c=1;
} else m |= 0x800000; } else m |= 0x800000;
val=m*ieee_table.e[c];
if(s) val = -val;
return val; val=m*ieee_table.e[c];
if(s) val = -val;
return val;
} }
unsigned long grib_ieee_nearest_smaller_to_long(double x) unsigned long grib_ieee_nearest_smaller_to_long(double x)
{ {
unsigned long l; unsigned long l;
unsigned long e; unsigned long e;
unsigned long m ; unsigned long m ;
unsigned long s; unsigned long s;
unsigned long mmin = 0x800000; unsigned long mmin = 0x800000;
double y,eps; double y,eps;
if(x == 0) return 0; if(x == 0) return 0;
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
l=grib_ieee_to_long(x); l=grib_ieee_to_long(x);
y=grib_long_to_ieee(l); y=grib_long_to_ieee(l);
if ( x < y ) { if ( x < y ) {
if ( x < 0 && -x < ieee_table.vmin ) { if ( x < 0 && -x < ieee_table.vmin ) {
l=0x80800000; l=0x80800000;
} else { } else {
e = (l & 0x7f800000) >> 23; e = (l & 0x7f800000) >> 23;
m = (l & 0x007fffff) | 0x800000; m = (l & 0x007fffff) | 0x800000;
s = l & 0x80000000; s = l & 0x80000000;
if ( m == mmin ) { if ( m == mmin ) {
/* printf("grib_ieee_nearest_smaller_to_long: m == mmin (0x%lX) e=%lu\n",m,e); */ /* printf("grib_ieee_nearest_smaller_to_long: m == mmin (0x%lX) e=%lu\n",m,e); */
e = s ? e : e-1; e = s ? e : e-1;
if (e<1) e=1; if (e<1) e=1;
if (e>254) e=254; if (e>254) e=254;
/* printf("grib_ieee_nearest_smaller_to_long: e=%lu \n",e); */ /* printf("grib_ieee_nearest_smaller_to_long: e=%lu \n",e); */
} }
eps=ieee_table.e[e]; eps=ieee_table.e[e];
/* printf("grib_ieee_nearest_smaller_to_long: x<y\n"); */ /* printf("grib_ieee_nearest_smaller_to_long: x<y\n"); */
l=grib_ieee_to_long(y-eps); l=grib_ieee_to_long(y-eps);
/* printf("grib_ieee_nearest_smaller_to_long: grib_ieee_to_long(y-eps)=0x%lX y=%.10e eps=%.10e x=%.10e\n",l,y,eps,x); */ /* printf("grib_ieee_nearest_smaller_to_long: grib_ieee_to_long(y-eps)=0x%lX y=%.10e eps=%.10e x=%.10e\n",l,y,eps,x); */
}
} else return l;
if (x<grib_long_to_ieee(l)) {
printf("grib_ieee_nearest_smaller_to_long: x=%.20e grib_long_to_ieee(0x%lX)=%.20e\n",x,l,grib_long_to_ieee(l));
Assert(x >= grib_long_to_ieee(l));
} }
} else return l;
if (x<grib_long_to_ieee(l)) { return l;
printf("grib_ieee_nearest_smaller_to_long: x=%.20e grib_long_to_ieee(0x%lX)=%.20e\n",x,l,grib_long_to_ieee(l));
Assert(x >= grib_long_to_ieee(l));
}
return l;
} }
int grib_nearest_smaller_ieee_float(double a,double* ret) { int grib_nearest_smaller_ieee_float(double a,double* ret) {
unsigned long l=0; unsigned long l=0;
if (!ieee_table.inited) init_ieee_table(); if (!ieee_table.inited) init_ieee_table();
if (a>ieee_table.vmax) return GRIB_INTERNAL_ERROR; if (a>ieee_table.vmax) return GRIB_INTERNAL_ERROR;
l=grib_ieee_nearest_smaller_to_long(a); l=grib_ieee_nearest_smaller_to_long(a);
*ret=grib_long_to_ieee(l); *ret=grib_long_to_ieee(l);
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
#else #else
@ -244,87 +244,87 @@ int grib_nearest_smaller_ieee_float(double a,double* ret) {
double grib_ieeefloat_error(double x) {return 0;} double grib_ieeefloat_error(double x) {return 0;}
double grib_long_to_ieee(unsigned long x) { double grib_long_to_ieee(unsigned long x) {
unsigned long s = x & 0x80000000; unsigned long s = x & 0x80000000;
unsigned long c = (x & 0x7f800000) >> 23; unsigned long c = (x & 0x7f800000) >> 23;
unsigned long m; unsigned long m;
double val; double val;
long e; long e;
if(x == 0) return 0; if(x == 0) return 0;
Assert(c != 255); Assert(c != 255);
if(c == 0) { m = x & 0x007fffff; e = -126 - 23; } if(c == 0) { m = x & 0x007fffff; e = -126 - 23; }
else { m = (x & 0x007fffff) | (1<<23); e = c - 127 - 23; } else { m = (x & 0x007fffff) | (1<<23); e = c - 127 - 23; }
val = m; val = m;
while(e < 0) { val /= 2.0; e++; } while(e < 0) { val /= 2.0; e++; }
while(e > 0) { val *= 2.0; e--; } while(e > 0) { val *= 2.0; e--; }
if(s) val = -val; if(s) val = -val;
return val; return val;
} }
int grib_nearest_smaller_ieee_float(double a,double* x) { int grib_nearest_smaller_ieee_float(double a,double* x) {
double e = grib_long_to_ieee(grib_ieee_to_long(a)); double e = grib_long_to_ieee(grib_ieee_to_long(a));
double b = a; double b = a;
/* printf("----> a=%g e=%g e-a=%g\n",a,e,e-a); */ /* printf("----> a=%g e=%g e-a=%g\n",a,e,e-a); */
if( e > b ) if( e > b )
{
unsigned long ub = grib_ieee_to_long(b);
unsigned long ue;
while(e>b)
{ {
/* printf("a=%g e=%g e-a=%g\n",a,e,e-a); */ unsigned long ub = grib_ieee_to_long(b);
a -= (e-a); unsigned long ue;
e = grib_long_to_ieee(grib_ieee_to_long(a)); while(e>b)
{
/* printf("a=%g e=%g e-a=%g\n",a,e,e-a); */
a -= (e-a);
e = grib_long_to_ieee(grib_ieee_to_long(a));
}
ue = grib_ieee_to_long(e);
Assert((ue-ub) == 1);
} }
ue = grib_ieee_to_long(e);
Assert((ue-ub) == 1);
}
Assert(b >= e); Assert(b >= e);
*x=e; *x=e;
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
unsigned long grib_ieee_to_long(double x) { unsigned long grib_ieee_to_long(double x) {
/* double y = x; */ /* double y = x; */
unsigned long s = 0; unsigned long s = 0;
unsigned long m; unsigned long m;
long p = 0; long p = 0;
unsigned long e = 0; unsigned long e = 0;
if(x == 0) return 0; if(x == 0) return 0;
if(x < 0) { s = 1; x = -x; } if(x < 0) { s = 1; x = -x; }
while(x < 2) { x *= 2; p--; } while(x < 2) { x *= 2; p--; }
while(x >= 2) { x /= 2; p++; } while(x >= 2) { x /= 2; p++; }
if(p > 127 ) { if(p > 127 ) {
/* Overflow */ /* Overflow */
e = 255; e = 255;
m = 0; m = 0;
} }
else if(p < -126) { else if(p < -126) {
/* int i; */ /* int i; */
e = 0; e = 0;
/* printf("p=%ld x=%g %ld\n",p,x,p+126+23); */ /* printf("p=%ld x=%g %ld\n",p,x,p+126+23); */
m = x * grib_power(p+126+23,2); m = x * grib_power(p+126+23,2);
} }
else { else {
e = p + 127; e = p + 127;
m = x * ( 1 << 23); m = x * ( 1 << 23);
m &= 0x007fffff; m &= 0x007fffff;
} }
m = (s << 31) | ( e << 23 ) | m; m = (s << 31) | ( e << 23 ) | m;
return m; return m;
} }
#endif #endif
@ -332,91 +332,91 @@ unsigned long grib_ieee_to_long(double x) {
#ifdef IEEE #ifdef IEEE
unsigned long grib_ieee64_to_long(double x) { unsigned long grib_ieee64_to_long(double x) {
unsigned long lval = 0; unsigned long lval = 0;
#if IEEE_LE #if IEEE_LE
unsigned char s[8]={0,}; unsigned char s[8]={0,};
unsigned char* buf=(unsigned char*)&x; unsigned char* buf=(unsigned char*)&x;
int j=0; int j=0;
for (j=7;j>=0;j--) for (j=7;j>=0;j--)
s[j]= *(buf++); s[j]= *(buf++);
memcpy(&lval,s,8); memcpy(&lval,s,8);
#elif IEEE_BE #elif IEEE_BE
memcpy(&lval,&x,8); memcpy(&lval,&x,8);
#endif #endif
return lval; return lval;
} }
double grib_long_to_ieee64(unsigned long x){ double grib_long_to_ieee64(unsigned long x){
double dval = 0.0; double dval = 0.0;
#if IEEE_LE #if IEEE_LE
unsigned char s[8]={0,}; unsigned char s[8]={0,};
unsigned char* buf=(unsigned char*)&x; unsigned char* buf=(unsigned char*)&x;
int j=0; int j=0;
for (j=7;j>=0;j--) for (j=7;j>=0;j--)
s[j]= *(buf++); s[j]= *(buf++);
memcpy(&dval,s,8); memcpy(&dval,s,8);
#elif IEEE_BE #elif IEEE_BE
memcpy(&dval,&x,8); memcpy(&dval,&x,8);
#else #else
Assert(!"Neither IEEE_LE nor IEEE_BE defined."); Assert(!"Neither IEEE_LE nor IEEE_BE defined.");
#endif #endif
return dval; return dval;
} }
int grib_ieee_decode_array(grib_context* c,unsigned char* buf,size_t nvals,int bytes,double* val) { int grib_ieee_decode_array(grib_context* c,unsigned char* buf,size_t nvals,int bytes,double* val) {
int err=0,i=0,j=0; int err=0,i=0,j=0;
unsigned char s[8]={0,}; unsigned char s[8]={0,};
float fval; float fval;
double* pval=val; double* pval=val;
switch (bytes) { switch (bytes) {
case 4: case 4:
for (i=0;i<nvals;i++) { for (i=0;i<nvals;i++) {
#if IEEE_LE #if IEEE_LE
for (j=3;j>=0;j--) for (j=3;j>=0;j--)
s[j]=*(buf++); s[j]=*(buf++);
memcpy(&fval,s,4); memcpy(&fval,s,4);
val[i]=(double)fval; val[i]=(double)fval;
#elif IEEE_BE #elif IEEE_BE
memcpy(&fval,buf,4); memcpy(&fval,buf,4);
val[i]=(double)fval; val[i]=(double)fval;
buf+=4; buf+=4;
#endif #endif
} }
break; break;
case 8: case 8:
for (i=0;i<nvals;i++) { for (i=0;i<nvals;i++) {
#if IEEE_LE #if IEEE_LE
for (j=7;j>=0;j--) for (j=7;j>=0;j--)
s[j]=*(buf++); s[j]=*(buf++);
memcpy(pval++,s,8); memcpy(pval++,s,8);
#elif IEEE_BE #elif IEEE_BE
memcpy(pval++,buf,8); memcpy(pval++,buf,8);
buf+=8; buf+=8;
#endif #endif
} }
break; break;
default: default:
grib_context_log(c,GRIB_LOG_ERROR, grib_context_log(c,GRIB_LOG_ERROR,
"grib_ieee_decode_array: %d bits not implemented",bytes*8); "grib_ieee_decode_array: %d bits not implemented",bytes*8);
return GRIB_NOT_IMPLEMENTED; return GRIB_NOT_IMPLEMENTED;
} }
return err; return err;
} }
#else #else
int grib_ieee_decode_array(grib_context* c,unsigned char* buf,size_t nvals,int bytes,double* val) { int grib_ieee_decode_array(grib_context* c,unsigned char* buf,size_t nvals,int bytes,double* val) {
int err=0,i=0; int err=0,i=0;
long bitr=0; long bitr=0;
for(i=0;i< nvals;i++) for(i=0;i< nvals;i++)
val[i] = grib_long_to_ieee(grib_decode_unsigned_long(buf,&bitr,bytes*8)); val[i] = grib_long_to_ieee(grib_decode_unsigned_long(buf,&bitr,bytes*8));
return err; return err;
} }
#endif #endif
@ -424,60 +424,60 @@ int grib_ieee_decode_array(grib_context* c,unsigned char* buf,size_t nvals,int b
#ifdef IEEE #ifdef IEEE
int grib_ieee_encode_array(grib_context* c,double* val,size_t nvals,int bytes, int grib_ieee_encode_array(grib_context* c,double* val,size_t nvals,int bytes,
unsigned char* buf) { unsigned char* buf) {
int err=0,i=0,j=0; int err=0,i=0,j=0;
unsigned char s4[4]; unsigned char s4[4];
unsigned char s8[8]; unsigned char s8[8];
float fval=0; float fval=0;
double *pval=val; double *pval=val;
switch (bytes) { switch (bytes) {
case 4: case 4:
for (i=0;i<nvals;i++) { for (i=0;i<nvals;i++) {
fval=(float)val[i]; fval=(float)val[i];
#if IEEE_LE #if IEEE_LE
memcpy(s4,&(fval),4); memcpy(s4,&(fval),4);
for (j=3;j>=0;j--) for (j=3;j>=0;j--)
*(buf++)=s4[j]; *(buf++)=s4[j];
#elif IEEE_BE #elif IEEE_BE
memcpy(buf,&(fval),4); memcpy(buf,&(fval),4);
buf+=4; buf+=4;
#endif #endif
} }
break; break;
case 8: case 8:
for (i=0;i<nvals;i++) { for (i=0;i<nvals;i++) {
#if IEEE_LE #if IEEE_LE
memcpy(s8,pval++,8); memcpy(s8,pval++,8);
for (j=7;j>=0;j--) for (j=7;j>=0;j--)
*(buf++)=s8[j]; *(buf++)=s8[j];
#elif IEEE_BE #elif IEEE_BE
memcpy(buf,pval++,8); memcpy(buf,pval++,8);
buf+=8; buf+=8;
#endif #endif
} }
break; break;
default: default:
grib_context_log(c,GRIB_LOG_ERROR, grib_context_log(c,GRIB_LOG_ERROR,
"grib_ieee_encode_array: %d bits not implemented",bytes*8); "grib_ieee_encode_array: %d bits not implemented",bytes*8);
return GRIB_NOT_IMPLEMENTED; return GRIB_NOT_IMPLEMENTED;
} }
return err; return err;
} }
#else #else
int grib_ieee_encode_array(grib_context* c,double* val,size_t nvals,int bytes, int grib_ieee_encode_array(grib_context* c,double* val,size_t nvals,int bytes,
unsigned char* buf) { unsigned char* buf) {
int err=0,i=0; int err=0,i=0;
long bitr=0; long bitr=0;
for(i=0;i< nvals;i++) for(i=0;i< nvals;i++)
grib_encode_unsigned_long(buf, grib_ieee_to_long(val[i]), &bitr, bytes*8); grib_encode_unsigned_long(buf, grib_ieee_to_long(val[i]), &bitr, bytes*8);
return err; return err;
} }
#endif #endif

92
src/grib_iterator_class_gaussian.c
View File

@ -87,42 +87,42 @@ static void init_class(grib_iterator_class* c)
static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j); static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j);
static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
grib_iterator_gaussian* self = (grib_iterator_gaussian*)i; grib_iterator_gaussian* self = (grib_iterator_gaussian*)i;
double *lats; double *lats;
double laf; /* latitude of first point in degrees */ double laf; /* latitude of first point in degrees */
double lal; /* latitude of last point in degrees */ double lal; /* latitude of last point in degrees */
long trunc; /* number of parallels between a pole and the equator */ long trunc; /* number of parallels between a pole and the equator */
long lai; long lai;
long jScansPositively=0; long jScansPositively=0;
int size; int size;
double start; double start;
unsigned long istart=0; unsigned long istart=0;
int ret = GRIB_SUCCESS; int ret = GRIB_SUCCESS;
const char* latofirst = grib_arguments_get_name(h,args,self->carg++); const char* latofirst = grib_arguments_get_name(h,args,self->carg++);
const char* latoflast = grib_arguments_get_name(h,args,self->carg++); const char* latoflast = grib_arguments_get_name(h,args,self->carg++);
const char* numtrunc = grib_arguments_get_name(h,args,self->carg++); const char* numtrunc = grib_arguments_get_name(h,args,self->carg++);
const char* s_jScansPositively = grib_arguments_get_name(h,args,self->carg++); const char* s_jScansPositively = grib_arguments_get_name(h,args,self->carg++);
if((ret = grib_get_double_internal(h,latofirst, &laf))) return ret; if((ret = grib_get_double_internal(h,latofirst, &laf))) return ret;
if((ret = grib_get_double_internal(h,latoflast, &lal))) return ret; if((ret = grib_get_double_internal(h,latoflast, &lal))) return ret;
if((ret = grib_get_long_internal(h,numtrunc,&trunc))) return ret; if((ret = grib_get_long_internal(h,numtrunc,&trunc))) return ret;
if((ret = grib_get_long_internal(h,s_jScansPositively,&jScansPositively))) if((ret = grib_get_long_internal(h,s_jScansPositively,&jScansPositively)))
return ret; return ret;
start=laf; start=laf;
size=trunc*2; size=trunc*2;
lats = grib_context_malloc(h->context,size*sizeof(double)); lats = grib_context_malloc(h->context,size*sizeof(double));
ret = grib_get_gaussian_latitudes(trunc, lats); ret = grib_get_gaussian_latitudes(trunc, lats);
if(ret != GRIB_SUCCESS) { if(ret != GRIB_SUCCESS) {
grib_context_log(h->context, GRIB_LOG_ERROR,"error %d calculating gaussian points",ret); grib_context_log(h->context, GRIB_LOG_ERROR,"error %d calculating gaussian points",ret);
return ret; return ret;
} }
/* /*
for(loi=(trunc*2)-1;loi>=0;loi--) for(loi=(trunc*2)-1;loi>=0;loi--)
if(fabs(lats[loi] - lal) < glatPrecision) break; if(fabs(lats[loi] - lal) < glatPrecision) break;
@ -131,24 +131,24 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
} }
*/ */
binary_search(lats,size-1,start,&istart); binary_search(lats,size-1,start,&istart);
Assert(istart < size); Assert(istart < size);
if (jScansPositively) { if (jScansPositively) {
for(lai=0;lai<self->nam;lai++) { for(lai=0;lai<self->nam;lai++) {
self->las[lai] = lats[istart--]; self->las[lai] = lats[istart--];
/*if (istart<0) istart=size-1; this condition is always FALSE -- 'istart' is unsigned long */ /*if (istart<0) istart=size-1; this condition is always FALSE -- 'istart' is unsigned long */
}
} else {
for(lai=0;lai<self->nam;lai++) {
self->las[lai] = lats[istart++];
if (istart>size-1) istart=0;
}
} }
} else {
for(lai=0;lai<self->nam;lai++) {
self->las[lai] = lats[istart++];
if (istart>size-1) istart=0;
}
}
grib_context_free(h->context,lats); grib_context_free(h->context,lats);
return ret; return ret;
} }
static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j) static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j)
@ -156,18 +156,18 @@ static void binary_search(double xx[], const unsigned long n, double x, unsigned
/*This routine works only on descending ordered arrays*/ /*This routine works only on descending ordered arrays*/
#define EPSILON 1e-3 #define EPSILON 1e-3
unsigned long ju,jm,jl; unsigned long ju,jm,jl;
jl=0; jl=0;
ju=n; ju=n;
while (ju-jl > 1) { while (ju-jl > 1) {
jm=(ju+jl) >> 1; jm=(ju+jl) >> 1;
if (fabs(x-xx[jm]) < EPSILON) { if (fabs(x-xx[jm]) < EPSILON) {
/* found something close enough. We're done */ /* found something close enough. We're done */
*j=jm; *j=jm;
return; return;
} }
if (x < xx[jm]) jl=jm; if (x < xx[jm]) jl=jm;
else ju=jm; else ju=jm;
} }
*j=jl; *j=jl;
} }

240
src/grib_nearest.c
View File

@ -21,84 +21,84 @@ int grib_nearest_find(
double inlat, double inlon, double inlat, double inlon,
unsigned long flags, unsigned long flags,
double* outlats,double* outlons, double* outlats,double* outlons,
double* values,double* distances, int* indexes, size_t *len) double* values, double* distances, int* indexes, size_t *len)
{ {
grib_nearest_class *c = nearest->cclass; grib_nearest_class *c = nearest->cclass;
Assert( flags <= (GRIB_NEAREST_SAME_GRID|GRIB_NEAREST_SAME_DATA|GRIB_NEAREST_SAME_POINT) ); Assert( flags <= (GRIB_NEAREST_SAME_GRID|GRIB_NEAREST_SAME_DATA|GRIB_NEAREST_SAME_POINT) );
while(c) while(c)
{ {
grib_nearest_class *s = c->super ? *(c->super) : NULL; grib_nearest_class *s = c->super ? *(c->super) : NULL;
if(c->find) { if(c->find) {
int ret=c->find(nearest,h,inlat,inlon,flags,outlats,outlons,values,distances,indexes,len); int ret=c->find(nearest,h,inlat,inlon,flags,outlats,outlons,values,distances,indexes,len);
if (ret!=GRIB_SUCCESS) { if (ret!=GRIB_SUCCESS) {
if (inlon>0) inlon-=360; if (inlon>0) inlon-=360;
else inlon+=360; else inlon+=360;
ret=c->find(nearest,h,inlat,inlon,flags,outlats,outlons,values,distances,indexes,len); ret=c->find(nearest,h,inlat,inlon,flags,outlats,outlons,values,distances,indexes,len);
} }
return ret; return ret;
} }
c = s; c = s;
} }
Assert(0); Assert(0);
return 0; return 0;
} }
/* For this one, ALL init are called */ /* For this one, ALL init are called */
static int init_nearest(grib_nearest_class* c,grib_nearest* i, grib_handle *h, grib_arguments* args) static int init_nearest(grib_nearest_class* c,grib_nearest* i, grib_handle *h, grib_arguments* args)
{ {
if(c) { if(c) {
int ret = GRIB_SUCCESS; int ret = GRIB_SUCCESS;
grib_nearest_class *s = c->super ? *(c->super) : NULL; grib_nearest_class *s = c->super ? *(c->super) : NULL;
if(!c->inited) if(!c->inited)
{ {
if(c->init_class) c->init_class(c); if(c->init_class) c->init_class(c);
c->inited = 1; c->inited = 1;
}
if(s) ret = init_nearest(s,i,h,args);
if(ret != GRIB_SUCCESS) return ret;
if(c->init) return c->init(i,h, args);
} }
if(s) ret = init_nearest(s,i,h,args); return GRIB_INTERNAL_ERROR;
if(ret != GRIB_SUCCESS) return ret;
if(c->init) return c->init(i,h, args);
}
return GRIB_INTERNAL_ERROR;
} }
int grib_nearest_init(grib_nearest* i, grib_handle *h, grib_arguments* args) int grib_nearest_init(grib_nearest* i, grib_handle *h, grib_arguments* args)
{ {
return init_nearest(i->cclass,i,h,args); return init_nearest(i->cclass,i,h,args);
} }
/* For this one, ALL destroy are called */ /* For this one, ALL destroy are called */
int grib_nearest_delete(grib_nearest *i) int grib_nearest_delete(grib_nearest *i)
{ {
grib_nearest_class *c = i->cclass; grib_nearest_class *c = i->cclass;
while(c) while(c)
{ {
grib_nearest_class *s = c->super ? *(c->super) : NULL; grib_nearest_class *s = c->super ? *(c->super) : NULL;
if(c->destroy) c->destroy(i); if(c->destroy) c->destroy(i);
c = s; c = s;
} }
return 0; return 0;
} }
void grib_binary_search(double xx[], const unsigned long n, double x, void grib_binary_search(double xx[], const unsigned long n, double x,
int *ju,int *jl) int *ju,int *jl)
{ {
size_t jm=0; size_t jm=0;
int ascending=0; int ascending=0;
*jl=0; *jl=0;
*ju=n; *ju=n;
ascending=(xx[n] >= xx[0]); ascending=(xx[n] >= xx[0]);
while (*ju-*jl > 1) { while (*ju-*jl > 1) {
jm=(*ju+*jl) >> 1; jm=(*ju+*jl) >> 1;
if ((x >= xx[jm]) == ascending) *jl=jm; if ((x >= xx[jm]) == ascending) *jl=jm;
else *ju=jm; else *ju=jm;
} }
} }
#define RADIAN(x) ((x) * acos(0.0) / 90.0) #define RADIAN(x) ((x) * acos(0.0) / 90.0)
double grib_nearest_distance(double radius,double lon1, double lat1, double lon2, double lat2){ double grib_nearest_distance(double radius,double lon1, double lat1, double lon2, double lat2){
@ -106,90 +106,90 @@ double grib_nearest_distance(double radius,double lon1, double lat1, double lon2
double rlat2=RADIAN(lat2); double rlat2=RADIAN(lat2);
double rlon1=lon1; double rlon1=lon1;
double rlon2=lon2; double rlon2=lon2;
double a; double a;
if (rlon1 >=360) rlon1-=360.0; if (rlon1 >=360) rlon1-=360.0;
rlon1=RADIAN(rlon1); rlon1=RADIAN(rlon1);
if (rlon2 >=360) rlon2-=360.0; if (rlon2 >=360) rlon2-=360.0;
rlon2=RADIAN(rlon2); rlon2=RADIAN(rlon2);
a = sin(rlat1)*sin(rlat2)+ cos(rlat1)*cos(rlat2)*cos(rlon2-rlon1); a = sin(rlat1)*sin(rlat2)+ cos(rlat1)*cos(rlat2)*cos(rlon2-rlon1);
if (a > 1 || a < -1 ) a=(int)a; if (a > 1 || a < -1 ) a=(int)a;
return radius*acos(a); return radius*acos(a);
} }
int grib_nearest_find_multiple(grib_handle* h,int is_lsm, int grib_nearest_find_multiple(grib_handle* h,int is_lsm,
double* inlats,double* inlons,long npoints, double* inlats,double* inlons,long npoints,
double* outlats,double* outlons, double* outlats,double* outlons,
double* values,double* distances, int* indexes) double* values,double* distances, int* indexes)
{ {
grib_nearest* nearest=0; grib_nearest* nearest=0;
double* pdistances=distances; double* pdistances=distances;
double* poutlats=outlats; double* poutlats=outlats;
double* poutlons=outlons; double* poutlons=outlons;
double* pvalues=values; double* pvalues=values;
int* pindexes=indexes; int* pindexes=indexes;
int idx=0,ii=0; int idx=0,ii=0;
double max,min; double max,min;
double qdistances[4]={0,}; double qdistances[4]={0,};
double qoutlats[4]={0,}; double qoutlats[4]={0,};
double qoutlons[4]={0,}; double qoutlons[4]={0,};
double qvalues[4]={0,}; double qvalues[4]={0,};
int qindexes[4]={0,}; int qindexes[4]={0,};
int ret=0; int ret=0;
long i=0; long i=0;
size_t len=4; size_t len=4;
int flags=GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA; int flags=GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA;
nearest=grib_nearest_new(h,&ret); nearest=grib_nearest_new(h,&ret);
if (ret!=GRIB_SUCCESS) return ret; if (ret!=GRIB_SUCCESS) return ret;
if (is_lsm) {
int noland=1;
for (i=0;i<npoints;i++) {
ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
qvalues,qdistances,qindexes,&len);
max=qdistances[0];
for (ii=0;ii<4;ii++) {
if (max<qdistances[ii]) {max=qdistances[ii];idx=ii;}
if (qvalues[ii] >= 0.5) noland=0;
}
min=max;
for (ii=0;ii<4;ii++) {
if ((min >= qdistances[ii]) && ( noland || (qvalues[ii] >= 0.5)) ) {
min = qdistances[ii];
idx=ii;
}
}
*poutlats=qoutlats[idx];poutlats++;
*poutlons=qoutlons[idx];poutlons++;
*pvalues=qvalues[idx];pvalues++;
*pdistances=qdistances[idx];pdistances++;
*pindexes=qindexes[idx];pindexes++;
if (is_lsm) {
int noland=1;
for (i=0;i<npoints;i++) {
ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
qvalues,qdistances,qindexes,&len);
max=qdistances[0];
for (ii=0;ii<4;ii++) {
if (max<qdistances[ii]) {max=qdistances[ii];idx=ii;}
if (qvalues[ii] >= 0.5) noland=0;
}
min=max;
for (ii=0;ii<4;ii++) {
if ((min >= qdistances[ii]) && ( noland || (qvalues[ii] >= 0.5)) ) {
min = qdistances[ii];
idx=ii;
} }
} } else {
*poutlats=qoutlats[idx];poutlats++; for (i=0;i<npoints;i++) {
*poutlons=qoutlons[idx];poutlons++; ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
*pvalues=qvalues[idx];pvalues++; qvalues,qdistances,qindexes,&len);
*pdistances=qdistances[idx];pdistances++; min=qdistances[0];
*pindexes=qindexes[idx];pindexes++; for (ii=0;ii<4;ii++) {
if ((min >= qdistances[ii])) {
} min = qdistances[ii];
} else { idx=ii;
for (i=0;i<npoints;i++) { }
ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons, }
qvalues,qdistances,qindexes,&len); *poutlats=qoutlats[idx];poutlats++;
min=qdistances[0]; *poutlons=qoutlons[idx];poutlons++;
for (ii=0;ii<4;ii++) { *pvalues=qvalues[idx];pvalues++;
if ((min >= qdistances[ii])) { *pdistances=qdistances[idx];pdistances++;
min = qdistances[ii]; *pindexes=qindexes[idx];pindexes++;
idx=ii;
} }
}
*poutlats=qoutlats[idx];poutlats++;
*poutlons=qoutlons[idx];poutlons++;
*pvalues=qvalues[idx];pvalues++;
*pdistances=qdistances[idx];pdistances++;
*pindexes=qindexes[idx];pindexes++;
} }
}
grib_nearest_delete(nearest); grib_nearest_delete(nearest);
return ret; return ret;
} }