00001
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "driver_utilities.h"
00033
00034
00039 void read_model_parameters( char* s_parm_name, int s_parm_relval) {
00040
00041 sprintf(msgStr,"Reading global parameters...");
00042 WriteMsg(msgStr,1);
00043 usrErr0(msgStr);
00044
00045 ReadGlobalParms(s_parm_name, s_parm_relval);
00046 sprintf(msgStr,"done.");
00047 WriteMsg(msgStr,1);
00048 usrErr(msgStr);
00049
00050 sprintf(msgStr,"Reading habitat-specific parameters...");
00051 WriteMsg(msgStr,1);
00052 usrErr0(msgStr);
00053
00054 ReadHabParms(s_parm_name, s_parm_relval);
00055 sprintf(msgStr,"done.");
00056 WriteMsg(msgStr,1);
00057 usrErr(msgStr);
00058
00059
00060 ReadModExperimParms(s_parm_name, s_parm_relval);
00061
00062 }
00063
00066 ViewParm* read_output_parms() {
00067
00068 ViewParm* v = readOutlist(Outlist_size);
00069 return v;
00070 }
00071
00075 ViewParm* readOutlist(int size)
00076 {
00077 int Npt;
00078 ViewParm* vp;
00079 FILE *cnfgFile;
00080
00081 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/RunParms/Model.outList",ModelPath,ProjName);
00082 else sprintf(modelFileName,"%s%s:RunParms:Model.outList",ModelPath,ProjName);
00083 cnfgFile = fopen(modelFileName, "r");
00084 if (cnfgFile == NULL) {
00085 sprintf(msgStr, "Error, can't open file %s", modelFileName);
00086 WriteMsg(msgStr, 1); exit(-1);}
00087
00088 Npt = readViewParms(cnfgFile,size, &vp);
00089 return vp;
00090
00091 }
00092
00093
00105 int readViewParms(FILE *vpFile, int size, ViewParm **ppview)
00106 {
00107 char test;
00108 char cmd;
00109 int index=0, ix, iy, iz, i, nArgs, precision, format=0, Ocmd;
00110 float x,y;
00111 ViewParm *vp, *view;
00112
00113 format = init_config_file(vpFile, '#', '*', '@', '~');
00114 view = (ViewParm *)nalloc(size * sizeof(ViewParm), "view");
00115 *ppview = view;
00116
00117 for (i = 0; i < size; i++) {
00118 vp = view+i;
00119 vp->scale.s = 1.0;
00120 vp->scale.o = 0.0;
00121 vp->step = 0;
00122 Ocmd=0;
00123 vp->nPoints = 0;
00124 vp->points = NULL;
00125 vp->maxPoints = 0;
00126 vp->mapType = 'N';
00127 precision = 1;
00128 setPrecision(vp,precision);
00129 }
00130
00131 while(1) {
00132 index = parse_packet(vpFile, &nArgs, &test);
00133 if( index == -1) break;
00134 if( index == -3) break;
00135 if(index >= size) {
00136 fprintf(stderr,"\n Read index Error in configuration: index out of range: %d ", index);
00137 break;
00138 }
00139 if( index == -2) break;
00140 else {
00141 vp = view+index;
00142 strncpy(vp->fileName,gCArg[0],23);
00143 vp->fileName[23]='\0';
00144 if( test == '*' ) setFlag(vp,ISARRAY,1);
00145 else if (test == '@') setFlag(vp,ISARRAY,0);
00146 else {
00147 sprintf(msgStr," %s: Syntax Error in configuration, test char = %c ",gCArg[0],test);
00148 usrErr(msgStr);
00149 break;
00150 }
00151 if(debug && gCArg[1][0] != 'U') {
00152 sprintf(msgStr,"\nReading Output Config for %s: Nargs = %d, CMD0 = %c index = %d",
00153 gCArg[0],nArgs,gCArg[1][0],index);
00154 WriteMsg(msgStr,1);
00155 }
00156 cmd = gCArg[1][0];
00157 switch (cmd) {
00158 case 'U':
00159 strncpy(vp->varUnits,gCArg[2],23);
00160 vp->varUnits[23]='\0';
00161 if(debug > 3) {
00162 sprintf(msgStr,"\n%s Output Config: Units= %s",gCArg[0],gCArg[2]);
00163 WriteMsg(msgStr,1);
00164 }
00165
00166 break;
00167 case 'O':
00168 if( isInteger(gCArg[2]) ) { ix = atoi(gCArg[2]); }
00169 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00170 if(debug) {
00171 sprintf(msgStr,"\n%s Output Config: O(%d)",gCArg[0],ix);
00172 WriteMsg(msgStr,1);
00173 }
00174 vp->step = ix; Ocmd=1;
00175
00176
00177 if (vp->step == CalMonOut && avg_Intvl != 0) {
00178 sprintf(msgStr,"\n***ERROR in Model.outList & Driver.parm integration: Driver.parm's avg_Intvl must be 0 for calendar-based Outstep of %s.",gCArg[0]);
00179 usrErr(msgStr);
00180 exit(-1);
00181 }
00182
00183 break;
00184 case 'A':
00185 if( isInteger(gCArg[2]) && nArgs > 2 ) { ix = atoi(gCArg[2]); }
00186 else {
00187 fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00188 WriteMsg(msgStr,1);
00189 }
00190 if(debug) {
00191 sprintf(msgStr,"\n%s Output Config: A()",gCArg[0]);
00192 WriteMsg(msgStr,1);
00193 }
00194 vp->mapType = 'H';
00195 if( Ocmd == 0 ) vp->step = 1;
00196 break;
00197 case 'G':
00198 if ( nArgs > 2 ) {
00199 if( isInteger(gCArg[2]) ) {
00200 vp->mapType = atoi(gCArg[2]);
00201 }
00202 else { vp->mapType =gCArg[2][0]; }
00203 }
00204 if ( nArgs > 3 )
00205 if( isInteger(gCArg[3]) ) {
00206 precision = atoi(gCArg[3]);
00207 }
00208 if ( nArgs > 4 ) strcpy(vp->fileName,gCArg[4]);
00209 if( precision > 1 && precision < 5 ) setPrecision(vp,precision);
00210 if( Ocmd == 0 ) vp->step = 1;
00211 break;
00212 case 'B':
00213 if( isFloat(gCArg[2]) && isFloat(gCArg[3]) && nArgs > 3 ) {
00214 vp->bounds.Vmax = atof(gCArg[2]);
00215 vp->bounds.Vmin = atof(gCArg[3]);
00216 }
00217 if(debug) {
00218 sprintf(msgStr,"\n Output Config: G(%.1f,%.1f)",
00219 vp->bounds.Vmax, vp->bounds.Vmin);
00220 WriteMsg(msgStr,1);
00221 }
00222 break;
00223 case 'P':
00224 if( isInteger(gCArg[2]) && isInteger(gCArg[3]) && nArgs > 3 ) {
00225 if( (i = (++(vp->nPoints)-1)) >= (vp->maxPoints)) make_more_points(vp,7);
00226 ix = atoi(gCArg[2]); iy = atoi(gCArg[3]);
00227 }
00228 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00229 if(debug) {
00230 sprintf(msgStr,"\n%s Output Config: P(%d,%d), index = %d",gCArg[0],ix,iy,i);
00231 WriteMsg(msgStr,1);
00232 }
00233 vp->points[i].x = ix;
00234 vp->points[i].y = iy;
00235 vp->points[i].type = 'p';
00236
00237
00238 if (vp->step == CalMonOut) {
00239 sprintf(msgStr,"\n***ERROR in Model.outList: sorry, but you can't have point-time series output with calendar-based Outstep of %s. Go yell at Carl!",gCArg[0]);
00240 usrErr(msgStr);
00241 exit(-1);
00242 }
00243 break;
00244 case 'W':
00245 if( isInteger(gCArg[2]) && isInteger(gCArg[3]) && isInteger(gCArg[4]) && nArgs > 4 ) {
00246 if( (i = (++(vp->nPoints)-1)) >= (vp->maxPoints)) make_more_points(vp,7);
00247 ix = atoi(gCArg[2]); iy = atoi(gCArg[3]); iz = atoi(gCArg[4]);
00248 }
00249 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00250 if( nArgs > 5 ) vp->points[i].format = gCArg[5][0];
00251 vp->points[i].type = 'w';
00252 vp->points[i].x = ix;
00253 vp->points[i].y = iy;
00254 vp->points[i].z = iz;
00255 if(debug) {
00256 sprintf(msgStr,"\n%s Output Config: W(%d,%d,%d,%c), index = %d",
00257 gCArg[0],ix,iy,iz,vp->points[i].format,i);
00258 WriteMsg(msgStr,1);
00259 }
00260 if( Ocmd == 0 ) vp->step = 1;
00261 break;
00262 case 'S': case 'C':
00263 if( gCArg[1][1] == 'm' || gCArg[1][0] == 'C' ) {
00264 if( isInteger(gCArg[2]) && isInteger(gCArg[3]) && isInteger(gCArg[4]) && nArgs > 5 ) {
00265 if( (i = (++(vp->nPoints)-1)) >= (vp->maxPoints)) make_more_points(vp,7);
00266 ix = atoi(gCArg[2]); iy = atoi(gCArg[3]); iz = atoi(gCArg[4]);
00267 vp->points[i].type = gCArg[5][0];
00268 }
00269 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00270 vp->points[i].x = ix;
00271 vp->points[i].y = iy;
00272 vp->points[i].z = iz;
00273 if(debug) {
00274 sprintf(msgStr,"\n%s Output Config: Sm(%d,%d,%d,%c), index = %d",
00275 gCArg[0],ix,iy,iz,gCArg[4][0],i);
00276 WriteMsg(msgStr,1);
00277 }
00278 if( Ocmd == 0 ) vp->step = 1;
00279 }
00280 else {
00281 if( nArgs > 2 ) {
00282 if( isFloat(gCArg[2]) ) { x = atof(gCArg[2]); }
00283 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00284 }
00285 if( nArgs > 3 ) {
00286 if( isFloat(gCArg[3]) ) { y = atof(gCArg[3]); }
00287 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00288 }
00289 if(debug) {
00290 sprintf(msgStr,"\n%s Output Config: S(%f,%f)",gCArg[0],x,y);
00291 WriteMsg(msgStr,1);
00292 }
00293 vp->scale.s = x;
00294 vp->scale.o = y;
00295 }
00296 break;
00297 }
00298 }
00299 }
00300 return size;
00301 }
00302
00303
00312 void make_more_points(ViewParm *vp, int ptIncr)
00313 {
00314 Point3D *new_pts;
00315 int i;
00316
00317 new_pts = (Point3D *) malloc( (vp->maxPoints+ptIncr) * sizeof(Point3D) );
00318 for(i=0; i < vp->maxPoints; i++) {
00319 new_pts[i].x = vp->points[i].x;
00320 new_pts[i].y = vp->points[i].y;
00321 new_pts[i].z = vp->points[i].z;
00322 new_pts[i].type = vp->points[i].type;
00323 }
00324 if( vp->points != NULL ) { free((char*)vp->points); }
00325 vp->points = new_pts;
00326 vp->maxPoints += ptIncr;
00327 if(debug) {
00328 sprintf(msgStr,"Made %d more points for %s for %d total points",ptIncr,vp->fileName,vp->maxPoints);
00329 WriteMsg(msgStr,1);
00330 }
00331 }
00332
00333
00344 void calc_maxmin(ViewParm *vp, void* Map, char Mtype, char* mName, int step)
00345 {
00346 int ix, iy;
00347 float vout[4], ftmp, fmax = -1000, fmin = 1000;
00348
00349 switch(Mtype) {
00350 case 'f' :
00351 for(ix=1; ix<=s0; ix++)
00352 for(iy=1; iy<=s1; iy++) {
00353 if( ON_MAP[T(ix,iy)] ) {
00354 if( (ftmp = ((float*)Map)[T(ix,iy)]) > fmax ) {
00355 fmax = ftmp;
00356 }
00357 if( ftmp < fmin ) fmin = ftmp;
00358 }
00359 }
00360 break;
00361 case 'i' : case 'd' :
00362 for(ix=1; ix<=s0; ix++)
00363 for(iy=1; iy<=s1; iy++) {
00364 if( ON_MAP[T(ix,iy)] ) {
00365 if((ftmp = ((int*)Map)[T(ix,iy)]) > fmax ) fmax = ftmp;
00366 if( ftmp < fmin ) fmin = ftmp;
00367 }
00368 }
00369 break;
00370 case 'c' :
00371 for(ix=1; ix<=s0; ix++)
00372 for(iy=1; iy<=s1; iy++) {
00373 if( ON_MAP[T(ix,iy)] ) {
00374 if( (ftmp = (float)((unsigned char*)Map)[T(ix,iy)]) > fmax ) fmax = ftmp;
00375 if( ftmp < fmin ) fmin = ftmp;
00376 }
00377 }
00378 break;
00379 }
00380 if(step==0) {
00381 vp->gScale.Vmax = fmax;
00382 vp->gScale.Vmin = fmin;
00383 } else {
00384 if( fmax > vp->gScale.Vmax ) vp->gScale.Vmax = fmax;
00385 if( fmin < vp->gScale.Vmin ) vp->gScale.Vmin = fmin;
00386 }
00387 vout[0] = fmax; vout[1] = vp->gScale.Vmax; vout[2] = fmin; vout[3] = vp->gScale.Vmin;
00388 Combine(vout,mName,4,kMAXMIN,step);
00389 }
00390
00392 void setup_grid() {
00393
00394 exgridsize(procnum,gbl_size,lcl_size,lcl_start);
00395
00396
00397 s0 = lcl_size[0];
00398 s1 = lcl_size[1];
00399
00400 gridSize = (s0+2)*(s1+2);
00401
00402 sprintf(msgStr,"\nGRID DATA::[ gsize: (%d, %d), lstart: (%d, %d), lend: (%d, %d), lsize: (%d, %d) ]\n",
00403 gbl_size[0], gbl_size[1], lcl_start[0],
00404 lcl_start[1], lcl_start[0]+lcl_size[0]-1,
00405 lcl_start[1]+lcl_size[1]-1, lcl_size[0], lcl_size[1] );
00406 WriteMsg(msgStr,1);
00407 sprintf(msgStr,"\nVP DATA:: size: (%d), Variable sizes: float: %d, int: %d, long: %d, double: %d\n",
00408 sizeof(ViewParm), sizeof(float),sizeof(int) ,sizeof(long) ,sizeof(double) );
00409 WriteMsg(msgStr,1);
00410
00411 gTempSize = gridSize*8;
00412 gTemp = (UCHAR*) nalloc(gTempSize,"gTemp");
00413
00414
00415
00416 floatOut = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"floatOut");
00417
00418
00419
00420 nc_dataRowRev = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"nc_dataRowRev");
00421 init_pvar(nc_dataRowRev,NULL,'f',0.0);
00422
00423
00424 }
00425
00426
00435 void write_output(int index, ViewParm* vp, void* Map,
00436 const char* filename, char Mtype, int step)
00437 {
00438 int i; Point3D point;
00439
00440
00441 if(vp->mapType != 'N' ) {
00442 write_map_file((char*)filename,Map,Mtype,index,vp->scale.s,vp->scale.o,getPrecision(vp),vp->mapType,vp->varUnits,vp->step,vp->fileName);
00443 if (debug > 3) calc_maxmin(vp,Map,Mtype,(char*)filename,step);
00444 }
00445 for(i=0; i< (vp->nPoints); i++ ) {
00446 point = vp->points[i];
00447 if (debug > 3) { sprintf(msgStr,"\nwriting Out: %s(%d): %c(), index = %d, step=%d\n", filename, i, point.type, index, step );
00448 WriteMsg(msgStr,1);
00449 }
00450
00451 switch( point.type ) {
00452 case 'm': case 'M': calc_maxmin( vp,Map,Mtype,(char*)filename,step); break;
00453 case 'w': quick_look( Map,(char*) filename, Mtype, point.x, point.y, point.z, point.format ); break;
00454 case 'a': case 's': case 'A': case 'S': print_loc_ave( &point, Map, Mtype, (char*)filename, step ); break;
00455 case 'p': print_point( vp, Map, Mtype, (char*)filename, step, i ); break;
00456 }
00457 }
00458 }
00459
00473 int write_map_file(const char* filename, VOIDP Map, char Mtype,
00474 int index, float scale_value, float offset_value,
00475 int bSize, unsigned char type, const char* thisvarUnits, int outstep, char* ncFileName)
00476 {
00477 int i, j, pathType=0;
00478 char ftype[10], gSize;
00479 float ftmp;
00480 SLONG itmp, imax, imin;
00481 UCHAR *mPtr;
00482 char s_mo[3], s_da[3];
00483 FILE* bFile;
00484
00485
00486
00487
00488
00489
00490
00491
00492 if (SimTime.mo[0]<10) sprintf(s_mo,"0%d", SimTime.mo[0]);
00493 else sprintf(s_mo,"%d", SimTime.mo[0]);
00494 if (SimTime.da[0]<10) sprintf(s_da,"0%d", SimTime.da[0]);
00495 else sprintf(s_da,"%d", SimTime.da[0]);
00496 sprintf(ftype,"%d%s%s\0",SimTime.yr[0],s_mo,s_da);
00497
00498 if( HDF4 && type=='H' ) sprintf(ftype,".hdf");
00499
00500 if( NCDF3 && type=='C' ) sprintf(ftype,".nc");
00501
00502 gSize = gridSize*bSize + 1;
00503 if( gSize > gTempSize ) {
00504 if(gTempSize) free((char*)gTemp);
00505 gTemp = (UCHAR*) nalloc( gSize, "gTemp" );
00506 gTempSize = gSize;
00507 }
00508
00509 if( H_OPSYS == UNIX ){
00510 if(check_for((char*)filename, "/", 0, CASE, END ) >= 0 ) {
00511 if(filename[0] != '/' ) pathType = 1;
00512 else pathType = 2;
00513 }
00514 }
00515 else { if(check_for ((char*)filename, ":", 0, CASE, END ) >= 0 ) { if( filename[0] == ':' ) pathType = 1; else pathType = 2; }}
00516
00517 if(type=='H') {
00518 switch(pathType) {
00519 case 0:
00520 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/HDF/%s%s\0",OutputPath,ProjName,filename,ftype);
00521 else sprintf(modelFileName,"%s%s:Output:HDF:%s%s\0",OutputPath,ProjName,filename,ftype);
00522 break;
00523 case 1:
00524 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/HDF/%s%s\0",OutputPath,ProjName,filename,ftype);
00525 else sprintf(modelFileName,"%s%s:Output:HDF:%s%s\0",OutputPath,ProjName,filename,ftype);
00526 break;
00527 case 2:
00528 sprintf(modelFileName,"%s.hdf\0",filename);
00529 break;
00530 }
00531 }
00532
00533 else if(type=='C') {
00534 switch(pathType) {
00535 case 0:
00536 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/NCDF/%s%s\0",OutputPath,ProjName,ncFileName,ftype);
00537 else sprintf(modelFileName,"%s%s:Output:NCDF:%s%s\0",OutputPath,ProjName,ncFileName,ftype);
00538 break;
00539 case 1:
00540 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/NCDF/%s%s\0",OutputPath,ProjName,ncFileName,ftype);
00541 else sprintf(modelFileName,"%s%s:Output:NCDF:%s%s\0",OutputPath,ProjName,ncFileName,ftype);
00542 break;
00543 case 2:
00544 sprintf(modelFileName,"%s.hdf\0",filename);
00545 break;
00546 }
00547 }
00548
00549 else if(type=='M') {
00550 switch(pathType) {
00551 case 0:
00552
00553 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/%s/%s%s\0",OutputPath,ProjName,filename,filename,ftype);
00554 else sprintf(modelFileName,"%s%s:Output:%s:%s%s\0",OutputPath,ProjName,filename,filename,ftype);
00555 break;
00556 case 1:
00557 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/%s/%s%s\0",OutputPath,ProjName,filename,filename,ftype);
00558 else sprintf(modelFileName,"%s%s:Output:%s:%s%s\0",OutputPath,ProjName,filename,filename,ftype);
00559 break;
00560 case 2:
00561 sprintf(modelFileName,"%s%s\0",filename,ftype);
00562 break;
00563 }
00564 }
00565 else {
00566 bSize = 1;
00567 switch(pathType) {
00568 case 0:
00569 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Animation%d/%s%s\0",OutputPath,ProjName,type,filename,ftype);
00570 else sprintf(modelFileName,"%s%s:Output:Animation%d:%s%s\0",OutputPath,ProjName,type,filename,ftype);
00571 break;
00572 case 1:
00573 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Animation%d/%s%s\0",OutputPath,ProjName,type,filename,ftype);
00574 else sprintf(modelFileName,"%s%s:Output:Animation%d:%s%s\0",OutputPath,ProjName,type,filename,ftype);
00575 break;
00576 case 2:
00577 sprintf(modelFileName,"%s%s\0",filename,ftype);
00578 break;
00579 }
00580 }
00581
00582 bSize = (bSize < 1) ? 1 : bSize;
00583 bSize = (bSize > 4) ? 4 : bSize;
00584 if( Mtype == 'c' ) bSize = 1;
00585
00586
00587
00588
00589
00590
00591 if (bSize == 4) {
00592
00593 for (i=0; i< (s0+2); i++) {
00594 for (j=0; j< (s1+2); j++) {
00595 if( ON_MAP[T(i,j)] ) {
00596 floatOut[T(i,j)] = (( (float*) Map)[T(i,j)]);
00597 }
00598 else {
00599 floatOut[T(i,j)] = offMap_float;
00600
00601 }
00602 }
00603 }
00604
00605 if ( debug >3) { sprintf(msgStr,"\n\nWriting Float Map %s: no scaling should be used (scale and offset values are ignored, but should be 1 & 0, respectively) scale_value = %f, offset_value = %f, index = %d, filename = %s, type = %c, bSize = %d, pathType = %d, type = %d, units = %s\n",
00606 filename,scale_value,offset_value,index,modelFileName,Mtype,bSize,pathType,type,thisvarUnits); WriteMsg(msgStr,1); }
00607
00608 writeFloatMap(modelFileName,filename,floatOut,bSize,type,index,thisvarUnits,outstep);
00609
00610 return(0);
00611 }
00612
00613
00614 if(bSize == 1) { imax = 255; imin = 0; }
00615 else if(bSize == 2) { imax = 255*128; imin = -imax; }
00616 else if(bSize == 3) { imax = 256*256*128; imin = -imax; }
00617 mPtr = gTemp;
00618
00619 if ( debug >3) { sprintf(msgStr,"\n\nWriting Map %s: scale_value = %f, offset_value = %f, index = %d, filename = %s, type = %c, bSize = %d, pathType = %d, type = %d\n",
00620 filename,scale_value,offset_value,index,modelFileName,Mtype,bSize,pathType,type); WriteMsg(msgStr,1); }
00621
00622 for (i=0; i<s0; i++) {
00623 for (j=0; j<s1; j++) {
00624 if( ON_MAP[T(i+1,j+1)] ) {
00625 switch(Mtype) {
00626 case 'l' :
00627 ftmp = ((( (double*) Map)[T(i+1,j+1)] - offset_value) / scale_value);
00628 itmp = (int) ftmp;
00629 itmp = (itmp > imax-1) ? imax-1 : itmp;
00630 itmp = (itmp < imin) ? imin : itmp;
00631 break;
00632 case 'f' :
00633 ftmp = ((( (float*) Map)[T(i+1,j+1)] - offset_value) / scale_value);
00634 itmp = (int) ftmp;
00635 itmp = (itmp > imax-1) ? imax-1 : itmp;
00636 itmp = (itmp < imin) ? imin : itmp;
00637 break;
00638 case 'd' : case 'i' :
00639 ftmp = ((( (int*) Map)[T(i+1,j+1)] - offset_value) / scale_value);
00640 itmp = (int) ftmp;
00641 itmp = (itmp > imax-1) ? imax-1 : itmp;
00642 itmp = (itmp < imin) ? imin : itmp;
00643 break;
00644 case 'c' :
00645 itmp = ((unsigned char*) Map)[T(i+1,j+1)];
00646 itmp = (itmp > imax-1) ? imax-1 : itmp;
00647 itmp = (itmp < imin) ? imin : itmp;
00648 break;
00649 }
00650 }
00651 else itmp = imax;
00652
00653
00654 enc_Nb(&mPtr,itmp,bSize);
00655
00656 }
00657 }
00658
00659 writeMap(modelFileName,gTemp,bSize,type,index);
00660 if(debug) quick_look(Map, (char*)filename, Mtype, dbgPt.x, dbgPt.y, 2,'E');
00661 return(0);
00662 }
00663
00664
00674 int read_map_file(const char *filename, VOIDP Map,
00675 unsigned char Mtype, float scale_value,
00676 float offset_value)
00677 {
00678 int i, j, k0, size;
00679 unsigned temp;
00680 UCHAR *tmpPtr, Dsize;
00681
00682 gridSize = (s0+2)*(s1+2);
00683 gTempSize = gridSize*8;
00684
00685
00686 size = gridSize*4 +1;
00687 if( size > gTempSize ) {
00688 if(gTemp) free((char*)gTemp);
00689 gTemp = (UCHAR*)nalloc( size, "gTemp" );
00690 gTempSize = size;
00691 }
00692
00693 if(debug>2) {
00694 sprintf(msgStr,"Reading %s\n",filename);
00695 usrErr(msgStr);
00696 WriteMsg(msgStr,1);
00697 }
00698
00699 Dsize = readMap(filename,gTemp);
00700
00701 if(debug) {
00702 sprintf(msgStr,"name = %s, proc = %d, scale_value = %f, offset_value = %f, size = %x\n ",
00703 filename,procnum,scale_value,offset_value,Dsize);
00704 WriteMsg(msgStr,1);
00705 }
00706 for (i=1; i<=s0; i++) {
00707 for (j=1; j<=s1; j++) {
00708 k0 = Dsize*((i-1)*lcl_size[1] + (j-1));
00709 tmpPtr = gTemp+k0;
00710 switch(Dsize) {
00711 case 1: temp = gTemp[k0]; break;
00712 case 2: temp = gTemp[k0]*256 + gTemp[k0+1]; break;
00713 case 3: temp = gTemp[k0]*65536 + gTemp[k0+1]*256 + gTemp[k0+2]; break;
00714 case 4: temp = gTemp[k0]*16777216 + gTemp[k0+1]*65536 + gTemp[k0+2]*256 + gTemp[k0+3]; break;
00715 default: fprintf(stderr,"ERROR, illegal size: %x\n",Dsize); temp = 0;
00716 }
00717 switch (Mtype) {
00718 case 'f' :
00719 ((float*)Map)[T(i,j)] = scale_value*((float)temp)+offset_value;
00720 break;
00721 case 'd' : case 'i' :
00722 ((int*)Map)[T(i,j)] = (int)(scale_value * (float)temp + offset_value);
00723 break;
00724 case 'c' :
00725 ((unsigned char*)Map)[T(i,j)] = (int)(scale_value * (unsigned char)temp);
00726 break;
00727 }
00728 }
00729 }
00730 if(debug) quick_look(Map, (char*)filename, Mtype, dbgPt.x, dbgPt.y, 2,'E');
00731 link_edges(Map,Mtype);
00732 fflush(stderr);
00733 return Dsize;
00734 }
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00761 void enc_Nb(UCHAR** sptr,SLONG value,int bytes)
00762 {
00763 UCHAR tmp[4];
00764 switch(bytes) {
00765 case 1:
00766 *(*sptr)++ = value;
00767 break;
00768 case 2:
00769 *(*sptr)++ = (value >> 8);
00770 *(*sptr)++ = value;
00771 break;
00772 case 3:
00773 tmp[0] = value;
00774 tmp[1] = (value >>= 8);
00775 tmp[2] = (value >> 8);
00776 *(*sptr)++ = tmp[2];
00777 *(*sptr)++ = tmp[1];
00778 *(*sptr)++ = tmp[0];
00779 break;
00780 case 4:
00781 tmp[0] = value;
00782 tmp[1] = (value >>= 8);
00783 tmp[2] = (value >>= 8);
00784 tmp[3] = (value >> 8);
00785 *(*sptr)++ = tmp[3];
00786 *(*sptr)++ = tmp[2];
00787 *(*sptr)++ = tmp[1];
00788 *(*sptr)++ = tmp[0];
00789 break;
00790 }
00791 }
00792
00802 void quick_look(void* Map, char* name, unsigned char Mtype,
00803 int xc, int yc, int range, unsigned char format)
00804 {
00805 int ymin, ymax, xmin, xmax, x, y, N0, N1;
00806 char* namePtr;
00807
00808 if (!on_this_proc(xc,yc)) return;
00809 x = xc - lcl_start[0] + 1;
00810 y = yc - lcl_start[1] + 1;
00811 xmin = x - range;
00812 xmax = x + range + 1;
00813 ymin = y - range;
00814 ymax = y + range + 1;
00815 xmax = iclip( xmax, 0, s0+2);
00816 xmin = iclip( xmin, 0, s0+2);
00817 ymax = iclip( ymax, 0, s1+2);
00818 ymin = iclip( ymin, 0, s1+2);
00819 N0 = xmax-xmin;
00820 N1 = ymax-ymin;
00821 namePtr = name_from_path(name);
00822 if (debug >3) {
00823 sprintf(msgStr,"Window to Map %s, s=%d, proc0 = %d, proc1 = %d; Corners(gbl): (%d,%d) (%d,%d)\n",
00824 namePtr,range,recpnum[0],recpnum[1],xc-range,yc-range,xc+range+1,yc+range+1);
00825 writeWindow(Map,namePtr,msgStr,xmin,ymin,N0,N1,istep,Mtype,format);
00826 }
00827
00828 }
00829
00830
00843 void writeWindow(void* fValue, char* vName, char* desc,
00844 int x0, int y0, int N0, int N1, int index,
00845 UCHAR Mtype, UCHAR format)
00846 {
00847 int mSize, size;
00848 char ctemp[200];
00849
00850 switch(Mtype) {
00851 case 'l': size = sizeof(double); break;
00852 case 'f': case 'L': case 'E': size = sizeof(float); break;
00853 case 'd': case 'i': size = sizeof(int); break;
00854 case 'c': size = sizeof(char); break;
00855 }
00856 if( (mSize=size*N0*N1) > gTempSize ) {
00857 if(gTemp != NULL) free((char*)gTemp);
00858 gTemp = (UCHAR*)malloc(gTempSize=mSize);
00859 }
00860
00861 Copy(((UCHAR*)fValue)+T(x0,y0)*size, gTemp, N1*size, N0, (s1+2)*size, N1*size);
00862 sprintf(ctemp,"(%d)WIND: %s",index,vName);
00863 writeSeries(gTemp,ctemp,desc,N0,N1,Mtype,format);
00864 }
00865
00872 int iclip ( int x0, int min, int max )
00873 {
00874 int rv;
00875 rv = ( x0 > max) ? max : x0;
00876 rv = ( rv < min) ? min : rv;
00877 return rv;
00878 }
00879
00880
00889 void Copy(void *src, void *dst, int w, int n, int sw, int dw)
00890 {
00891 int i;
00892 for(i=0; i<n; i++)
00893 memcpy( ((UCHAR*)dst)+i*dw, ((UCHAR*)src)+i*sw, w );
00894 }
00895
00896
00897
00898
00900 void set_env_vars(void) {
00901
00902 FILE *infile;
00903 char filename[100], ch;
00904 static long start;
00905 int i, maxLen =200;
00906
00907 ModelPath = getenv("ModelPath");
00908 ProjName = getenv("ProjName");
00909 DriverPath = getenv("DriverPath");
00910 OS_TYPE = getenv("OSTYPE");
00911
00912 sprintf(msgStr,"OS type = %s ",OS_TYPE);
00913 usrErr(msgStr);
00914
00915 sprintf(msgStr,"Project Name = %s ",ProjName);
00916 usrErr(msgStr);
00917
00918 }
00919
00920
00929 void print_point(ViewParm *vp, void* Map, char Mtype,
00930 char* mName, int tIndex, int vpindex)
00931 {
00932 int x, y;
00933 Point3D *pt;
00934
00935 pt = vp->points + vpindex;
00936 x = pt->x - lcl_start[0];
00937 y = pt->y - lcl_start[1];
00938 if( (x >= 0) && (x < s0) && (y >= 0) && (y < s1) ) {
00939 if(tIndex==0) { pt->z = PListIndex; pSeries[PListIndex].Length=0;
00940 }
00941 else PListIndex = pt->z;
00942 if(tIndex==0) {
00943
00944 if(PListIndex >= MAX_PTSERIES ) fatal("Too many Point series.");
00945 else if( pSeries[PListIndex].data == NULL ) pSeries[PListIndex].data = (float*) malloc( (N_iter*sizeof(float)/vp->step)+2 );
00946 if( pSeries[PListIndex].data == NULL ) usrErr("out of Memory for Spatial Timeseries.");
00947 pSeries[PListIndex].Loc.x = pt->x;
00948 pSeries[PListIndex].Loc.y = pt->y;
00949 pSeries[PListIndex].outstep = vp->step;
00950 strcpy(pSeries[PListIndex].name,mName);
00951 if(debug >2) {
00952 sprintf(msgStr,"\nSetup Pointlist %d for %s(%d), step=%d, point=(%d,%d)\n", PListIndex, mName, vpindex, vp->step, x, y );
00953 WriteMsg(msgStr,1);
00954 }
00955 }
00956 switch(Mtype) {
00957 case 'f' : pSeries[PListIndex].data[ pSeries[PListIndex].Length++ ] = ((float*)Map)[T(x+1,y+1)]; break;
00958 case 'i' : pSeries[PListIndex].data[ pSeries[PListIndex].Length++ ] = ((int*)Map)[T(x+1,y+1)]; break;
00959 case 'c' : pSeries[PListIndex].data[ pSeries[PListIndex].Length++ ] = ((unsigned char*)Map)[T(x+1,y+1)];
00960 }
00961 if(debug >3) {
00962 sprintf(msgStr,"\nWriting Point %d for %s(%d), point=(%d,%d), value = %f, index = %d\n",
00963 PListIndex, mName, vpindex, x, y, ((float*)Map)[T(x+1,y+1)], pSeries[PListIndex].Length-1 );
00964 WriteMsg(msgStr,1);
00965 }
00966 PListIndex++;
00967 }
00968 }
00969
00975 void read_map_dims(const char* filename)
00976 {
00977 FILE *file;
00978
00979 if (debug>3) { sprintf(msgStr,"Getting map dims: %s",filename); usrErr(msgStr);}
00980 if(Lprocnum == 1) {
00981 sprintf(modelFileName,"%s/%s/Data/Map_head/%s",ModelPath,ProjName,filename);
00982 file = fopen(modelFileName,"r");
00983 if(file==NULL) {
00984 fprintf(stderr,"Unable to open map header file %s.\n",modelFileName);
00985 fflush(stderr);
00986 Exit(0);
00987 }
00988 scan_forward(file,"cols:");
00989 fscanf(file,"%d",&gbl_size[1]);
00990 sprintf(msgStr,"cols = %d\n",gbl_size[1]); WriteMsg(msgStr,1);
00991 scan_forward(file,"rows:");
00992 fscanf(file,"%d",&gbl_size[0]);
00993 sprintf(msgStr,"rows = %d\n",gbl_size[0]); WriteMsg(msgStr,1);
00994 fclose(file);
00995 }
00996 broadcastInt(gbl_size);
00997 broadcastInt(gbl_size+1);
00998 s0 = gbl_size[0];
00999 s1 = gbl_size[1];
01000 }
01001
01008 void init_pvar(VOIDP Map, UCHAR* mask, unsigned char Mtype,float iv)
01009 {
01010 int i0, i1;
01011
01012 switch(Mtype) {
01013 case 'b' :
01014 for(i0=0; i0<=numBasn; i0++) {
01015 ((double*)Map)[i0] = iv;
01016 }
01017 break;
01018 case 'l' :
01019 for(i0=0; i0<=s0+1; i0++)
01020 for(i1=0; i1<=s1+1; i1++) {
01021 if(mask==NULL) ((double*)Map)[T(i0,i1)] = iv;
01022 else if ( mask[T(i0,i1)] == 0 ) ((double*)Map)[T(i0,i1)] = 0;
01023 else ((double*)Map)[T(i0,i1)] = iv;
01024 }
01025 break;
01026 case 'f' :
01027 for(i0=0; i0<=s0+1; i0++)
01028 for(i1=0; i1<=s1+1; i1++) {
01029 if(mask==NULL) ((float*)Map)[T(i0,i1)] = iv;
01030 else if ( mask[T(i0,i1)] == 0 ) ((float*)Map)[T(i0,i1)] = 0;
01031 else ((float*)Map)[T(i0,i1)] = iv;
01032 }
01033 break;
01034 case 'i' : case 'd' :
01035 for(i0=0; i0<=s0+1; i0++)
01036 for(i1=0; i1<=s1+1; i1++) {
01037 if(mask==NULL) ((int*)Map)[T(i0,i1)] = (int)iv;
01038 else if ( mask[T(i0,i1)] == 0 ) ((int*)Map)[T(i0,i1)] = 0;
01039 else ((int*)Map)[T(i0,i1)] = (int)iv;
01040 }
01041 break;
01042 case 'c' :
01043 for(i0=0; i0<=s0+1; i0++)
01044 for(i1=0; i1<=s1+1; i1++) {
01045 if(mask==NULL) ((unsigned char*)Map)[T(i0,i1)] = (unsigned char) '\0' + (int) iv;
01046 else if ( mask[T(i0,i1)] == 0 ) ((unsigned char*)Map)[T(i0,i1)] = '\0';
01047 else ((unsigned char*)Map)[T(i0,i1)] = (unsigned char) '\0' + (int) iv;
01048 }
01049 break;
01050
01051 default :
01052 printf(" in default case\n");
01053 break;
01054 }
01055 }
01056
01065 void print_loc_ave(Point3D *vt, void* Map, char Mtype, char* mName, int tIndex)
01066 {
01067 int ix, iy, x, y, x0, y0, range, xmax, ymax, xmin, ymin, type = 11;
01068 float vout[4];
01069 double sum = 0, normCoef = 0;
01070
01071 if(Mtype != 'f') { fprintf(stderr,"Warning: Wrong type in Loc Ave: %s(%c)\n",mName,Mtype); return; }
01072 x0 = vt->x;
01073 y0 = vt->y;
01074 range = vt->z;
01075 x = x0 - lcl_start[0];
01076 y = y0 - lcl_start[1];
01077 xmin = x - range;
01078 xmax = x + range;
01079 ymin = y - range;
01080 ymax = y + range;
01081 xmax = (xmax > (s0)) ? (s0) : xmax;
01082 xmin = (xmin < 0 ) ? 0 : xmin;
01083 ymax = (ymax > (s1)) ? (s1) : ymax;
01084 ymin = (ymin < 0 ) ? 0 : ymin;
01085
01086 for(ix=xmin; ix<xmax; ix++) {
01087 for(iy=ymin; iy<ymax; iy++) {
01088 if( ON_MAP[T(ix+1,iy+1)] ) sum += ((float*)Map)[T(ix+1,iy+1)];
01089 normCoef += 1.0;
01090 }
01091 }
01092
01093 vout[0] = sum; vout[1] = normCoef;
01094
01095 switch (vt->type) {
01096 case 's' : Combine(vout,mName,1,kSUM,tIndex); break;
01097 case 'a' : Combine(vout,mName,2,kAVE,tIndex); break;
01098 case 'A' : Combine(vout,mName,2,kAVECUM,tIndex); break;
01099 case 'S' : Combine(vout,mName,1,kSUMCUM,tIndex); break;
01100 }
01101 }
01102
01103
01105 void PTS_SetFields (PTSeries* thisStruct, int ix, int iy,
01106 int index, int format, int col)
01107 {
01108 float ptstep = 1.0, *value;
01109 int n0;
01110
01111 thisStruct->fX = ix;
01112 thisStruct->fY = iy;
01113 if(thisStruct->fValue) free((char*)thisStruct->fValue);
01114
01115 thisStruct->fValue = readSeriesCol(modelFileName,format,index,&n0,&ptstep,col);
01116
01117 gNPtTs = thisStruct->fNpt = n0;
01118 gTS = thisStruct->fTS = ptstep;
01119 value = thisStruct->fValue;
01120 if ( debug>2 ) { sprintf(msgStr,"\nPTS_SetFields: TSVals: %f %f %f, TS: %f, NPt: %d \n",value[0],value[1],value[2],ptstep,n0); WriteMsg(msgStr,1); }
01121 }
01122
01124 void PTS_Free(PTSeries* thisStruct)
01125 {
01126 if(thisStruct->fValue) free((char*)thisStruct->fValue);
01127 }
01128
01129
01131 void PTS_CopyFields (PTSeries* thisStruct, PTSeries pV)
01132 {
01133 thisStruct->fX = pV.fX;
01134 thisStruct->fY = pV.fY;
01135 thisStruct->fValue = pV.fValue;
01136 }
01137
01139 PTSeriesList* PTSL_Init(int nSlots, int format)
01140 {
01141 int i;
01142 PTSeriesList* thisStruct;
01143 thisStruct = (PTSeriesList*) malloc(sizeof(PTSeriesList));
01144 thisStruct->fList = (PTSeries*) malloc(sizeof(PTSeries)*nSlots);
01145 for(i=0; i<nSlots; i++) thisStruct->fList[i].fValue = NULL;
01146 thisStruct->fNSlots = nSlots;
01147 thisStruct->fNSlotsUsed = 0;
01148 thisStruct->fNptTS = 0;
01149 thisStruct->fFormat = format;
01150 return(thisStruct);
01151 }
01152
01154 void PTSL_AddpTSeries(PTSeriesList* thisStruct, int x, int y, int index, int seriesNum, int col )
01155 {
01156 PTSeries* valListtemp;
01157 int i, newSlots;
01158
01159 if( seriesNum >= thisStruct->fNSlots ) {
01160 newSlots = Max(20,(int)(thisStruct->fNSlots * 0.2));
01161 thisStruct->fNSlots += newSlots;
01162 valListtemp = (PTSeries*) malloc(sizeof(PTSeries)*(thisStruct->fNSlots));
01163 for(i=0; i<thisStruct->fNSlotsUsed; i++) {
01164 valListtemp[i].fNpt = thisStruct->fNptTS;
01165 PTS_CopyFields(valListtemp+i,thisStruct->fList[i]);
01166 }
01167 for(i=thisStruct->fNSlotsUsed; i<thisStruct->fNSlots; i++) valListtemp[i].fValue = NULL;
01168 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01169 thisStruct->fList = valListtemp;
01170 }
01171
01172 PTS_SetFields( thisStruct->fList + seriesNum, x, y, index, thisStruct->fFormat, col );
01173 thisStruct->fNSlotsUsed = seriesNum;
01174 }
01175
01177 void PTSL_Free(PTSeriesList* thisStruct)
01178 {
01179 int i;
01180 if(thisStruct == NULL) return;
01181 for(i=0; i<thisStruct->fNSlotsUsed; i++) PTS_Free( thisStruct->fList +i );
01182 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01183 free((char*)thisStruct);
01184 thisStruct = NULL;
01185 }
01186
01188 float PTSL_GetInterpolatedValue0(PTSeriesList* thisStruct, int x, int y, int step)
01189 {
01190 int i, dx, dy, my_debug;
01191 PTSeries pV;
01192 float wpow;
01193 float weight, InterpValue=0.0, distance=0.0;
01194 int r;
01195
01196 wpow = GP_IDW_pow;
01197
01198
01199 for(i=0; i<thisStruct->fNSlotsUsed; i++) {
01200 pV = thisStruct->fList[i];
01201 dx = (pV.fX-(x+lcl_start[0]));
01202 dy = (pV.fY-(y+lcl_start[1]));
01203 r = (dx*dx + dy*dy);
01204 if( r == 0 ) return pV.fValue[step];
01205 weight = (wpow == 1.0) ? (1.0)/r : 1.0/(r*r);
01206
01207
01208
01209 InterpValue += pV.fValue[step]*weight;
01210 distance += weight;
01211 }
01212 if (distance > 0) InterpValue /= distance;
01213 else InterpValue = (thisStruct->fList[0]).fValue[step];
01214 return (float) InterpValue;
01215 }
01216
01218 void PTSL_ReadLists(PTSeriesList* thisStruct, const char *ptsFileName,
01219 int index, float* timeStep, int* nPtTS, int col )
01220 {
01221 char pathname[150], infilename[60], ss[201], ret = '\n';
01222 FILE* cfile;
01223 int ix, iy, tst, sCnt=0;
01224 if( Lprocnum == 1 ) {
01225 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Data/%s.pts",ModelPath,ProjName,ptsFileName);
01226 else sprintf(modelFileName,"%s%s:Data:%s.pts",ModelPath,ProjName,ptsFileName);
01227 cfile = fopen(modelFileName,"r");
01228 if(cfile==NULL) {fprintf(stdout,"\nERROR: Unable to open timeseries file %s\n",modelFileName); Exit(1); }
01229 else { sprintf(msgStr,"\nReading file %s\n",modelFileName); WriteMsg(msgStr,1); }
01230 if (debug > 2) {sprintf(msgStr,"Reading %s timeseries, column %d",modelFileName, col); usrErr(msgStr);}
01231
01232
01233
01234 fgets(ss,200,cfile);
01235 }
01236 while(1) {
01237 if( Lprocnum == 1 ) {
01238 tst = fscanf(cfile,"%d",&ix);
01239 if( tst > 0 ) {
01240 fscanf(cfile,"%d",&iy);
01241 tst = fscanf(cfile,"%s",infilename);
01242 sprintf(modelFileName,"%s/%s/Data/%s",ModelPath,ProjName,infilename);
01243 }
01244 }
01245 broadcastInt(&tst);
01246 if(tst<1) break;
01247 broadcastInt(&ix);
01248 broadcastInt(&iy);
01249 PTSL_AddpTSeries(thisStruct, ix, iy, index, sCnt, col);
01250 sCnt++;
01251 }
01252 if( Lprocnum == 1 ) fclose(cfile);
01253 *nPtTS = gNPtTs;
01254 *timeStep = gTS;
01255 }
01256
01258 void PTSL_CreatePointMap(PTSeriesList* pList,void* Map, unsigned char Mtype,
01259 int step, float scale)
01260 {
01261 int i0, i1;
01262
01263 switch(Mtype) {
01264 case 'f' :
01265 for(i0=0; i0<=s0+1; i0++)
01266 for(i1=0; i1<=s1+1; i1++)
01267 ((float*)Map)[T(i0,i1)] = (ON_MAP[T(i0,i1)]) ? PTSL_GetInterpolatedValue0(pList,i0,i1,step)*scale : 0.0 ;
01268 break;
01269 case 'i' : case 'd' :
01270 for(i0=0; i0<=s0+1; i0++)
01271 for(i1=0; i1<=s1+1; i1++)
01272 ((int*)Map)[T(i0,i1)] = (int) ( (ON_MAP[T(i0,i1)]) ? PTSL_GetInterpolatedValue0(pList,i0,i1,step)*scale : 0 );
01273 break;
01274 case 'c' :
01275 for(i0=0; i0<=s0+1; i0++)
01276 for(i1=0; i1<=s1+1; i1++)
01277 ((unsigned char*)Map)[T(i0,i1)] = (unsigned char) '\0' + (int) ( (ON_MAP[T(i0,i1)]) ? PTSL_GetInterpolatedValue0(pList,i0,i1,step)*scale : 0 );
01278 break;
01279 }
01280 }
01281
01283 float PTSL_GetInterpolatedValue1(PTSeriesList* thisStruct, int x, int y, int step)
01284 {
01285 int i, dx, dy;
01286 PTSeries pV;
01287 RPointList *pList;
01288 float wpow;
01289 double weight, InterpValue=0.0, distance=0.0;
01290 long r;
01291
01292 pList = RPL_Init( 20 );
01293 wpow = GP_IDW_pow;
01294 for(i=0; i<thisStruct->fNSlotsUsed; i++) {
01295 pV = thisStruct->fList[i];
01296 dx = (pV.fX-x);
01297 dy = (pV.fY-y);
01298 r = dx*dx + dy*dy;
01299 if( r == 0 ) return pV.fValue[step];
01300 RPL_AddrPoint(pList, dx, dy, r, pV.fValue[step]);
01301 weight = (wpow == 1.0) ? ((double)1.0)/r : pow(r,-wpow);
01302 InterpValue += pV.fValue[step]*weight;
01303 distance += weight;
01304 }
01305 RPL_Sort(pList);
01306 RPL_Free(pList);
01307 return (float) InterpValue;
01308 }
01309
01311 void RP_SetFields (RPoint* thisStruct, int ix, int iy, float r, float value)
01312 {
01313 thisStruct->fX = ix;
01314 thisStruct->fY = iy;
01315 thisStruct->fr = r;
01316 thisStruct->fValue = value;
01317 }
01318
01320 void RP_CopyFields (RPoint* thisStruct, RPoint* that)
01321 {
01322 thisStruct->fX = that->fX;
01323 thisStruct->fY = that->fY;
01324 thisStruct->fr = that->fr;
01325 thisStruct->fValue = that->fValue;
01326 }
01327
01329 void RP_SwapFields (RPoint* thisStruct, RPoint* that)
01330 {
01331 RPoint temp;
01332 RP_CopyFields (&temp, thisStruct);
01333 RP_CopyFields (thisStruct, that);
01334 RP_CopyFields (that, &temp);
01335 }
01336
01338 RPointList * RPL_Init(int nPoints)
01339 {
01340 RPointList *thisStruct;
01341 thisStruct = (RPointList*) malloc(sizeof(RPointList));
01342 thisStruct->fList = (RPoint*) malloc(sizeof(RPoint)*nPoints);
01343 thisStruct->fNPt = nPoints;
01344 thisStruct->fNPtUsed = 0;
01345 return thisStruct;
01346 }
01347
01349 void RPL_AddrPoint(RPointList *thisStruct, int x, int y, float r, float value)
01350 {
01351 RPoint* pointListtemp;
01352 int i, newPoints;
01353
01354 if( thisStruct->fNPtUsed >= thisStruct->fNPt ) {
01355 newPoints = 20;
01356 thisStruct->fNPt += newPoints;
01357 pointListtemp = (RPoint*) malloc(sizeof(RPoint)*(thisStruct->fNPt));
01358 for(i=0; i<thisStruct->fNPtUsed; i++) {
01359 RP_CopyFields(pointListtemp+i,thisStruct->fList+i);
01360 }
01361 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01362 thisStruct->fList = pointListtemp;
01363 }
01364
01365 RP_SetFields( thisStruct->fList + thisStruct->fNPtUsed, x, y, r, value);
01366 thisStruct->fNPtUsed++;
01367 }
01368
01370 void RPL_Sort(RPointList *thisStruct)
01371 {
01372 int i, go=1;
01373
01374 while(go) {
01375 go = 0;
01376 for(i=1; i<thisStruct->fNPtUsed; i++) if(thisStruct->fList[i-1].fr > thisStruct->fList[i].fr) {
01377 RP_SwapFields (thisStruct->fList+i-1, thisStruct->fList+i);
01378 go = 1;
01379 }
01380 }
01381 }
01382
01383
01385 void RPL_Free(RPointList* thisStruct)
01386 {
01387 if(thisStruct == NULL) return;
01388 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01389 free((char*)thisStruct);
01390 thisStruct = NULL;
01391 }
01392
01393
01405 void setFlag(ViewParm* vp, UINT mask, int value )
01406 {
01407 if(value) vp->flags |= mask; else vp->flags &= ~mask;
01408 }
01409
01411 int getFlag(ViewParm* vp, UINT mask)
01412 {
01413 if( vp->flags & mask ) return 1; else return 0;
01414 }
01415
01417 void setPrecision(ViewParm* vp, int value )
01418 {
01419 if(value/3) vp->flags |= PMASK2; else vp->flags &= ~PMASK2;
01420 if((value-1)%2) vp->flags |= PMASK1; else vp->flags &= ~PMASK1;
01421 }
01422
01424 int getPrecision(ViewParm* vp)
01425 {
01426 int rv = 1;
01427 if(vp->flags & PMASK2) rv +=2;
01428 if(vp->flags & PMASK1) rv +=1;
01429 return rv;
01430 }
01431
01443 int check_for(char *s0, const char *s1, int start, int cs, int rp)
01444 {
01445
01446
01447
01448
01449
01450 int i, j=0, k=-1, Len1 = strlen(s1), Len0 = strlen(s0);
01451 char t1, t2;
01452
01453 while(k<0) {
01454 k=0;
01455 for(i=start; i< Len0; ++i) {
01456 t1 = s0[i];
01457 t2 = s1[j];
01458 if(cs==NOCASE) { t1 = tolower(t1); t2 = tolower(t2); }
01459 if (t1 == t2) j++;
01460 else {j=0; k=0;}
01461 if(j==Len1) { k=1; break; }
01462 }
01463 }
01464 if(k<=0) return(-1);
01465 else if(rp==BEG) return(i-Len1+1);
01466 else return(i+1);
01467 }
01468
01469
01475 int isInteger (char *target_str) {
01476
01477 int i=-1,first_num=0;
01478 char ch;
01479
01480 while( (ch=target_str[++i]) != '\0' ) {
01481 if( isdigit(ch) ) first_num=1;
01482 if( (ch=='-' || ch=='+') && first_num ) return(0);
01483 if( !( isspace(ch) || isdigit(ch) || ch=='-' || ch=='+') ) return(0);
01484 }
01485 return(1);
01486 }
01487
01493 int isFloat (char *target_str) {
01494
01495 int i=-1,first_num=0;
01496 char ch;
01497
01498 while( (ch=target_str[++i]) != '\0' ) {
01499 if( isdigit(ch) ) first_num=1;
01500 if( (ch=='-' || ch=='+') && first_num ) return(0);
01501 if( !( isspace(ch) || isdigit(ch) || ch=='-' || ch=='+' || ch=='.' || toupper(ch) == 'E') ) return(0);
01502 }
01503 return(1);
01504 }
01505
01506
01519 void WriteMsg(const char* msg, int wh) {
01520 wh = 1;
01521 if(Driver_outfile) fprintf(Driver_outfile,"%s\n",msg);
01522 else fprintf(stdout,"%s\n",msg);
01523 fflush(stdout);
01524 }
01525
01530 void usrErr0(const char* dString)
01531 {
01532 fprintf(stderr,"%s", dString);
01533 fflush(stderr);
01534 }
01535
01540 void usrErr(const char* dString)
01541 {
01542 fprintf(stderr,"%s\n", dString);
01543 fflush(stderr);
01544 }
01545
01546
01548 void Exit(int code) { exit(code); }
01549
01551 void fatal(const char *msg)
01552 {
01553 printf("%s",msg);
01554 fflush(stdout);
01555 Exit(9);
01556 }
01557
01558
01559
01560
01561
01583 double julday(int mon, int day, int year, int h, int mi, double se)
01584 {
01585 long m = mon, d = day, y = year;
01586 long c, ya, j;
01587 double seconds = h * 3600.0 + mi * 60 + se;
01588
01589 if (m > 2)
01590 m -= 3;
01591 else {
01592 m += 9;
01593 --y;
01594 }
01595 c = y / 100L;
01596 ya = y - (100L * c);
01597 j = (146097L * c) / 4L + (1461L * ya) / 4L + (153L * m + 2L) / 5L + d + 1721119L;
01598 if (seconds < 12 * 3600.0) {
01599 j--;
01600 seconds += 12.0 * 3600.0;
01601 }
01602 else {
01603 seconds = seconds - 12.0 * 3600.0;
01604 }
01605 return (j + (seconds / 3600.0) / 24.0);
01606 }
01607
01630 void calcdate(double jd, int *m, int *d, int *y, int *h, int *mi, double *sec)
01631 {
01632 static int ret[4];
01633
01634 long j = (long)jd;
01635 double tmp, frac = jd - j;
01636
01637 if (frac >= 0.5) {
01638 frac = frac - 0.5;
01639 j++;
01640 }
01641 else {
01642 frac = frac + 0.5;
01643 }
01644
01645 ret[3] = (j + 1L) % 7L;
01646 j -= 1721119L;
01647 *y = (4L * j - 1L) / 146097L;
01648 j = 4L * j - 1L - 146097L * *y;
01649 *d = j / 4L;
01650 j = (4L * *d + 3L) / 1461L;
01651 *d = 4L * *d + 3L - 1461L * j;
01652 *d = (*d + 4L) / 4L;
01653 *m = (5L * *d - 3L) / 153L;
01654 *d = 5L * *d - 3 - 153L * *m;
01655 *d = (*d + 5L) / 5L;
01656 *y = 100L * *y + j;
01657 if (*m < 10)
01658 *m += 3;
01659 else {
01660 *m -= 9;
01661 *y += 1;
01662 }
01663 tmp = 3600.0 * (frac * 24.0);
01664 *h = (int) (tmp / 3600.0);
01665 tmp = tmp - *h * 3600.0;
01666
01667 *mi = (int) (tmp / 60.0);
01668 *sec = tmp - *mi * 60.0;
01669 }
01670
01675 float FMOD(float x, float y) {
01676 return (float)fmod((double)x, (double)y);
01677 }
01678
01683 void getFloat(FILE* inFile, const char* lString, float* fValPtr)
01684 {
01685 int test, fSize;
01686 UCHAR iEx=0;
01687 if(lString) scan_forward(inFile,lString);
01688 test = fscanf(inFile,"%f",fValPtr);
01689 if(test != 1) {
01690 fprintf(stderr,"Read Error (%d):",test);
01691 if(lString) fprintf(stderr,"%s\n",lString);
01692 *fValPtr = 0.0;
01693 iEx = 1;
01694 }
01695 if (iEx) exit(0);
01696 }
01697
01702 void getChar(FILE* inFile, const char* lString, char* cValPtr)
01703 {
01704 int test;
01705 UCHAR iEx=0;
01706 if(lString) scan_forward(inFile,lString);
01707 test = fscanf(inFile,"%c",cValPtr);
01708 if(test != 1) {
01709 fprintf(stderr,"Read Error (%d):",test);
01710 if(lString) fprintf(stderr,"%s\n",lString);
01711 iEx = 1;
01712 *cValPtr = '\0';
01713 }
01714 if (iEx) exit(0);
01715 }
01716
01717
01722 void getString(FILE* inFile, const char *lString, char *inString)
01723 {
01724 int test;
01725 UCHAR iEx=0;
01726 if(lString) scan_forward(inFile,lString);
01727 test = fscanf(inFile,"%s",inString);
01728 if(test != 1) {
01729 fprintf(stderr,"Read Error (%d):",test);
01730 if(lString) fprintf(stderr,"%s\n",lString);
01731 iEx = 1;
01732 }
01733 if (iEx) exit(0);
01734 }
01735
01740 void getInt(FILE* inFile, const char* lString, int* iValPtr)
01741 {
01742 int test;
01743 UCHAR iEx=0;
01744 if(lString)
01745 scan_forward(inFile,lString);
01746 test = fscanf(inFile,"%d",iValPtr);
01747 if(test != 1) {
01748 fprintf(stderr,"Read Error (%d):",test);
01749 if(lString) fprintf(stderr,"%s\n",lString);
01750 *iValPtr = 0;
01751 iEx = 1;
01752 }
01753 if (iEx) exit(0);
01754 }
01755
01760 int find_char(FILE *infile, char tchar)
01761 {
01762 char test = '1', cchar = '#';
01763 int in_c=0;
01764
01765 while( test != EOF ) {
01766 test = fgetc(infile);
01767 if(test == tchar && in_c==0 ) { return 1; }
01768 if(test == cchar ) in_c=1;
01769 if(test == '\n' ) in_c=0;
01770 }
01771 return(0);
01772 }
01773
01777 int Round(float x)
01778 {
01779 int i = (int)x;
01780 return (i);
01781 }
01782
01788 char *Scip(char *s, char SYM )
01789 {
01790 if(*s == SYM ) return ++s;
01791 while(*s != SYM && *s != EOS ) s++;
01792 if(*s != EOS) return ++s;
01793 else return s;
01794 }
01795
01796
01800 int skip_white(FILE *infile)
01801 {
01802 int ch;
01803
01804 while( isspace(ch=fgetc(infile)) ) {;}
01805 if(ch==EOF) return 0;
01806 ungetc(ch,infile);
01807 return 1;
01808 }
01809
01814 int scan_forward ( FILE *infile, const char *tstring)
01815 {
01816 int sLen, i, cnt=0;
01817 char Input_string[100], test;
01818
01819 sLen = strlen(tstring);
01820 while( ( test = fgetc(infile) ) != EOF ) {
01821 for(i=0; i<(sLen-1); i++)
01822 Input_string[i] = Input_string[i+1];
01823 Input_string[sLen-1] = test;
01824 Input_string[sLen] = '\0';
01825 if(++cnt >= sLen) {
01826 test = strcmp(Input_string,tstring);
01827 if( test == 0 ) return 1;
01828 }
01829 }
01830 return(-1);
01831 }
01832
01833
01837 char* name_from_path(char* name)
01838 {
01839 char* namePtr; int i, slen;
01840 char dirCh;
01841
01842 namePtr = name;
01843 dirCh = '/';
01844 slen = strlen(name);
01845
01846 for(i=0; i<slen; i++) {
01847 if( name[slen-i-1] == dirCh ) { namePtr = name+slen-i; break; }
01848 }
01849 return namePtr;
01850 }
01851
01852
01857 VOIDP nalloc(unsigned mem_size, const char var_name[])
01858 {
01859 VOIDP rp;
01860
01861
01862 if(mem_size == 0) return(NULL);
01863 rp = (VOIDP)malloc( mem_size );
01864 total_memory += mem_size;
01865 fasync(stderr);
01866 if( rp == NULL ) {
01867 fprintf(stderr,"Sorry, out of memory(%d): %s\n",mem_size,var_name);
01868 Exit(0);
01869 }
01870 fmulti(stderr);
01871 return(rp);
01872 }
01873
01875 float Normal(float mean,float sd)
01876 {
01877 int table_loc ;
01878 double rand_num ;
01879 float sign ;
01880 float interval = .1 ;
01881 int n_interval = 40 ;
01882 float high, low ;
01883 float return_val = 0.0;
01884
01885 sign = SMDRAND(0.0, 1.0) ;
01886 sign = (sign < .5 ? -1 : 1) ;
01887 rand_num = SMDRAND(.5, 1.0);
01888 low = gRTable[0] ;
01889 for (table_loc=1; table_loc<n_interval; table_loc++)
01890 {
01891 high = gRTable[table_loc] ;
01892 if (high > rand_num + .5)
01893 {
01894 return_val = mean + sd * (sign * interval * (table_loc - 1 +
01895 (rand_num+.5-low)/(high-low))) ;
01896 return(return_val) ;
01897 }
01898 low = high ;
01899 }
01900 return(return_val) ;
01901 }
01902
01904 int Poisson(float mu)
01905 {
01906 int ix;
01907 float f0, r, Lp = 1, xf = 1, p=0;
01908
01909 p = f0 = Exp(-mu);
01910 r = SMDRAND(0.0,1.0);
01911
01912 for( ix=0; ix<500; ix++) {
01913 if( r < p ) return ix;
01914 Lp *= mu;
01915 if(ix>0) xf *= ix;
01916 p += (Lp*f0)/xf;
01917 }
01918 return ix;
01919 }
01920
01921
01922
01924 void link_edges(VOIDP Map, unsigned char Mtype)
01925 {
01926 switch(Mtype) {
01927 case 'f' :
01928 exchange_borders((byte*)Map,sizeof(float)); break;
01929 case 'd' : case 'i' :
01930 exchange_borders((byte*)Map,sizeof(int)); break;
01931 case 'c' :
01932 exchange_borders((byte*)Map,sizeof(UCHAR)); break;
01933 }
01934 }
01935
01937 void setup_platform() {
01938
01939 exparam(&env);
01940 procnum = env.procnum;
01941 exgridsplit(env.nprocs, 2, nprocs);
01942 if( exgridinit(2, nprocs) < 0) { usrErr("Failed to Setup Grid"); exit(0); }
01943 exgridcoord(procnum, recpnum);
01944 }
01945
01946
01947
01948