00001
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "driver_utilities.h"
00029
00030
00035 void read_model_parameters( char* s_parm_name, int s_parm_relval) {
00036
00037 sprintf(msgStr,"Reading global parameters...");
00038 WriteMsg(msgStr,1);
00039 usrErr0(msgStr);
00040
00041 ReadGlobalParms(s_parm_name, s_parm_relval);
00042 sprintf(msgStr,"done.");
00043 WriteMsg(msgStr,1);
00044 usrErr(msgStr);
00045
00046 sprintf(msgStr,"Reading habitat-specific parameters...");
00047 WriteMsg(msgStr,1);
00048 usrErr0(msgStr);
00049
00050 ReadHabParms(s_parm_name, s_parm_relval);
00051 sprintf(msgStr,"done.");
00052 WriteMsg(msgStr,1);
00053 usrErr(msgStr);
00054
00055 }
00056
00059 ViewParm* read_output_parms() {
00060
00061 ViewParm* v = readOutlist(Outlist_size);
00062 return v;
00063 }
00064
00068 ViewParm* readOutlist(int size)
00069 {
00070 int Npt;
00071 ViewParm* vp;
00072 FILE *cnfgFile;
00073
00074 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/RunParms/Model.outList",ModelPath,ProjName);
00075 else sprintf(modelFileName,"%s%s:RunParms:Model.outList",ModelPath,ProjName);
00076 cnfgFile = fopen(modelFileName, "r");
00077 if (cnfgFile == NULL) {
00078 sprintf(msgStr, "Error, can't open file %s", modelFileName);
00079 WriteMsg(msgStr, 1); exit(-1);}
00080
00081 Npt = readViewParms(cnfgFile,size, &vp);
00082 return vp;
00083
00084 }
00085
00086
00098 int readViewParms(FILE *vpFile, int size, ViewParm **ppview)
00099 {
00100 char test;
00101 char cmd;
00102 int index=0, ix, iy, iz, i, nArgs, precision, format=0, Ocmd;
00103 float x,y;
00104 ViewParm *vp, *view;
00105
00106 format = init_config_file(vpFile, '#', '*', '@', '~');
00107 view = (ViewParm *)nalloc(size * sizeof(ViewParm), "view");
00108 *ppview = view;
00109
00110 for (i = 0; i < size; i++) {
00111 vp = view+i;
00112 vp->scale.s = 1.0;
00113 vp->scale.o = 0.0;
00114 vp->step = 0;
00115 Ocmd=0;
00116 vp->nPoints = 0;
00117 vp->points = NULL;
00118 vp->maxPoints = 0;
00119 vp->mapType = 'N';
00120 precision = 1;
00121 setPrecision(vp,precision);
00122 }
00123
00124 while(1) {
00125 index = parse_packet(vpFile, &nArgs, &test);
00126 if( index == -1) break;
00127 if( index == -3) break;
00128 if(index >= size) {
00129 fprintf(stderr,"\n Read index Error in configuration: index out of range: %d ", index);
00130 break;
00131 }
00132 if( index == -2) break;
00133 else {
00134 vp = view+index;
00135 strncpy(vp->fileName,gCArg[0],23);
00136 vp->fileName[23]='\0';
00137 if( test == '*' ) setFlag(vp,ISARRAY,1);
00138 else if (test == '@') setFlag(vp,ISARRAY,0);
00139 else {
00140 sprintf(msgStr," %s: Syntax Error in configuration, test char = %c ",gCArg[0],test);
00141 usrErr(msgStr);
00142 break;
00143 }
00144 if(debug) {
00145 sprintf(msgStr,"\nReading Output Config for %s: Nargs = %d, CMD0 = %c index = %d",
00146 gCArg[0],nArgs,gCArg[1][0],index);
00147 WriteMsg(msgStr,1);
00148 }
00149 cmd = gCArg[1][0];
00150 switch (cmd) {
00151 case 'O':
00152 if( isInteger(gCArg[2]) ) { ix = atoi(gCArg[2]); }
00153 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00154 if(debug) {
00155 sprintf(msgStr,"\n%s Output Config: O(%d)",gCArg[0],ix);
00156 WriteMsg(msgStr,1);
00157 }
00158 vp->step = ix; Ocmd=1;
00159
00160
00161 if (vp->step == CalMonOut && avg_Intvl != 0) {
00162 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]);
00163 usrErr(msgStr);
00164 exit(-1);
00165 }
00166
00167 break;
00168 case 'A':
00169 if( isInteger(gCArg[2]) && nArgs > 2 ) { ix = atoi(gCArg[2]); }
00170 else {
00171 fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00172 WriteMsg(msgStr,1);
00173 }
00174 if(debug) {
00175 sprintf(msgStr,"\n%s Output Config: A()",gCArg[0]);
00176 WriteMsg(msgStr,1);
00177 }
00178 vp->mapType = 'H';
00179 if( Ocmd == 0 ) vp->step = 1;
00180 break;
00181 case 'G':
00182 if ( nArgs > 2 ) {
00183 if( isInteger(gCArg[2]) ) {
00184 vp->mapType = atoi(gCArg[2]);
00185 }
00186 else { vp->mapType =gCArg[2][0]; }
00187 }
00188 if ( nArgs > 3 )
00189 if( isInteger(gCArg[3]) ) {
00190 precision = atoi(gCArg[3]);
00191 }
00192 if ( nArgs > 4 ) strcpy(vp->fileName,gCArg[4]);
00193 if( precision > 1 && precision < 5 ) setPrecision(vp,precision);
00194 if( Ocmd == 0 ) vp->step = 1;
00195 break;
00196 case 'B':
00197 if( isFloat(gCArg[2]) && isFloat(gCArg[3]) && nArgs > 3 ) {
00198 vp->bounds.Vmax = atof(gCArg[2]);
00199 vp->bounds.Vmin = atof(gCArg[3]);
00200 }
00201 if(debug) {
00202 sprintf(msgStr,"\n Output Config: G(%.1f,%.1f)",
00203 vp->bounds.Vmax, vp->bounds.Vmin);
00204 WriteMsg(msgStr,1);
00205 }
00206 break;
00207 case 'P':
00208 if( isInteger(gCArg[2]) && isInteger(gCArg[3]) && nArgs > 3 ) {
00209 if( (i = (++(vp->nPoints)-1)) >= (vp->maxPoints)) make_more_points(vp,7);
00210 ix = atoi(gCArg[2]); iy = atoi(gCArg[3]);
00211 }
00212 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00213 if(debug) {
00214 sprintf(msgStr,"\n%s Output Config: P(%d,%d), index = %d",gCArg[0],ix,iy,i);
00215 WriteMsg(msgStr,1);
00216 }
00217 vp->points[i].x = ix;
00218 vp->points[i].y = iy;
00219 vp->points[i].type = 'p';
00220
00221
00222 if (vp->step == CalMonOut) {
00223 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]);
00224 usrErr(msgStr);
00225 exit(-1);
00226 }
00227 break;
00228 case 'W':
00229 if( isInteger(gCArg[2]) && isInteger(gCArg[3]) && isInteger(gCArg[4]) && nArgs > 4 ) {
00230 if( (i = (++(vp->nPoints)-1)) >= (vp->maxPoints)) make_more_points(vp,7);
00231 ix = atoi(gCArg[2]); iy = atoi(gCArg[3]); iz = atoi(gCArg[4]);
00232 }
00233 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00234 if( nArgs > 5 ) vp->points[i].format = gCArg[5][0];
00235 vp->points[i].type = 'w';
00236 vp->points[i].x = ix;
00237 vp->points[i].y = iy;
00238 vp->points[i].z = iz;
00239 if(debug) {
00240 sprintf(msgStr,"\n%s Output Config: W(%d,%d,%d,%c), index = %d",
00241 gCArg[0],ix,iy,iz,vp->points[i].format,i);
00242 WriteMsg(msgStr,1);
00243 }
00244 if( Ocmd == 0 ) vp->step = 1;
00245 break;
00246 case 'S': case 'C':
00247 if( gCArg[1][1] == 'm' || gCArg[1][0] == 'C' ) {
00248 if( isInteger(gCArg[2]) && isInteger(gCArg[3]) && isInteger(gCArg[4]) && nArgs > 5 ) {
00249 if( (i = (++(vp->nPoints)-1)) >= (vp->maxPoints)) make_more_points(vp,7);
00250 ix = atoi(gCArg[2]); iy = atoi(gCArg[3]); iz = atoi(gCArg[4]);
00251 vp->points[i].type = gCArg[5][0];
00252 }
00253 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00254 vp->points[i].x = ix;
00255 vp->points[i].y = iy;
00256 vp->points[i].z = iz;
00257 if(debug) {
00258 sprintf(msgStr,"\n%s Output Config: Sm(%d,%d,%d,%c), index = %d",
00259 gCArg[0],ix,iy,iz,gCArg[4][0],i);
00260 WriteMsg(msgStr,1);
00261 }
00262 if( Ocmd == 0 ) vp->step = 1;
00263 }
00264 else {
00265 if( nArgs > 2 ) {
00266 if( isFloat(gCArg[2]) ) { x = atof(gCArg[2]); }
00267 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00268 }
00269 if( nArgs > 3 ) {
00270 if( isFloat(gCArg[3]) ) { y = atof(gCArg[3]); }
00271 else fprintf(stderr," %s: Syntax Error in configuration ",gCArg[0]);
00272 }
00273 if(debug) {
00274 sprintf(msgStr,"\n%s Output Config: S(%f,%f)",gCArg[0],x,y);
00275 WriteMsg(msgStr,1);
00276 }
00277 vp->scale.s = x;
00278 vp->scale.o = y;
00279 }
00280 break;
00281 }
00282 }
00283 }
00284 return size;
00285 }
00286
00287
00296 void make_more_points(ViewParm *vp, int ptIncr)
00297 {
00298 Point3D *new_pts;
00299 int i;
00300
00301 new_pts = (Point3D *) malloc( (vp->maxPoints+ptIncr) * sizeof(Point3D) );
00302 for(i=0; i < vp->maxPoints; i++) {
00303 new_pts[i].x = vp->points[i].x;
00304 new_pts[i].y = vp->points[i].y;
00305 new_pts[i].z = vp->points[i].z;
00306 new_pts[i].type = vp->points[i].type;
00307 }
00308 if( vp->points != NULL ) { free((char*)vp->points); }
00309 vp->points = new_pts;
00310 vp->maxPoints += ptIncr;
00311 if(debug) {
00312 sprintf(msgStr,"Made %d more points for %s for %d total points",ptIncr,vp->fileName,vp->maxPoints);
00313 WriteMsg(msgStr,1);
00314 }
00315 }
00316
00317
00328 void calc_maxmin(ViewParm *vp, void* Map, char Mtype, char* mName, int step)
00329 {
00330 int ix, iy;
00331 float vout[4], ftmp, fmax = -1000, fmin = 1000;
00332
00333 switch(Mtype) {
00334 case 'f' :
00335 for(ix=1; ix<=s0; ix++)
00336 for(iy=1; iy<=s1; iy++) {
00337 if( ON_MAP[T(ix,iy)] ) {
00338 if( (ftmp = ((float*)Map)[T(ix,iy)]) > fmax ) {
00339 fmax = ftmp;
00340 }
00341 if( ftmp < fmin ) fmin = ftmp;
00342 }
00343 }
00344 break;
00345 case 'i' : case 'd' :
00346 for(ix=1; ix<=s0; ix++)
00347 for(iy=1; iy<=s1; iy++) {
00348 if( ON_MAP[T(ix,iy)] ) {
00349 if((ftmp = ((int*)Map)[T(ix,iy)]) > fmax ) fmax = ftmp;
00350 if( ftmp < fmin ) fmin = ftmp;
00351 }
00352 }
00353 break;
00354 case 'c' :
00355 for(ix=1; ix<=s0; ix++)
00356 for(iy=1; iy<=s1; iy++) {
00357 if( ON_MAP[T(ix,iy)] ) {
00358 if( (ftmp = (float)((unsigned char*)Map)[T(ix,iy)]) > fmax ) fmax = ftmp;
00359 if( ftmp < fmin ) fmin = ftmp;
00360 }
00361 }
00362 break;
00363 }
00364 if(step==0) {
00365 vp->gScale.Vmax = fmax;
00366 vp->gScale.Vmin = fmin;
00367 } else {
00368 if( fmax > vp->gScale.Vmax ) vp->gScale.Vmax = fmax;
00369 if( fmin < vp->gScale.Vmin ) vp->gScale.Vmin = fmin;
00370 }
00371 vout[0] = fmax; vout[1] = vp->gScale.Vmax; vout[2] = fmin; vout[3] = vp->gScale.Vmin;
00372 Combine(vout,mName,4,kMAXMIN,step);
00373 }
00374
00376 void setup_grid() {
00377
00378 exgridsize(procnum,gbl_size,lcl_size,lcl_start);
00379
00380
00381 s0 = lcl_size[0];
00382 s1 = lcl_size[1];
00383
00384 gridSize = (s0+2)*(s1+2);
00385
00386 sprintf(msgStr,"\nGRID DATA::[ gsize: (%d, %d), lstart: (%d, %d), lend: (%d, %d), lsize: (%d, %d) ]\n",
00387 gbl_size[0], gbl_size[1], lcl_start[0],
00388 lcl_start[1], lcl_start[0]+lcl_size[0]-1,
00389 lcl_start[1]+lcl_size[1]-1, lcl_size[0], lcl_size[1] );
00390 WriteMsg(msgStr,1);
00391 sprintf(msgStr,"\nVP DATA:: size: (%d), Variable sizes: float: %d, int: %d, long: %d, double: %d\n",
00392 sizeof(ViewParm), sizeof(float),sizeof(int) ,sizeof(long) ,sizeof(double) );
00393 WriteMsg(msgStr,1);
00394
00395 gTempSize = gridSize*8;
00396 gTemp = (UCHAR*) nalloc(gTempSize,"gTemp");
00397 }
00398
00399
00408 void write_output(int index, ViewParm* vp, void* Map,
00409 const char* filename, char Mtype, int step)
00410 {
00411 int i; Point3D point;
00412
00413
00414 if(vp->mapType != 'N' ) {
00415 write_map_file((char*)filename,Map,Mtype,index,vp->scale.s,vp->scale.o,getPrecision(vp),vp->mapType);
00416 if (debug > 3) calc_maxmin(vp,Map,Mtype,(char*)filename,step);
00417 }
00418 for(i=0; i< (vp->nPoints); i++ ) {
00419 point = vp->points[i];
00420 if (debug > 3) { sprintf(msgStr,"\nwriting Out: %s(%d): %c(), index = %d, step=%d\n", filename, i, point.type, index, step );
00421 WriteMsg(msgStr,1);
00422 }
00423
00424 switch( point.type ) {
00425 case 'm': case 'M': calc_maxmin( vp,Map,Mtype,(char*)filename,step); break;
00426 case 'w': quick_look( Map,(char*) filename, Mtype, point.x, point.y, point.z, point.format ); break;
00427 case 'a': case 's': case 'A': case 'S': print_loc_ave( &point, Map, Mtype, (char*)filename, step ); break;
00428 case 'p': print_point( vp, Map, Mtype, (char*)filename, step, i ); break;
00429 }
00430 }
00431 }
00432
00443 int write_map_file(const char* filename, VOIDP Map, char Mtype,
00444 int index, float scale_value, float offset_value,
00445 int bSize, unsigned char type)
00446 {
00447 int i, j, pathType=0;
00448 char ftype[10], gSize;
00449 float ftmp;
00450 SLONG itmp, imax, imin;
00451 UCHAR *mPtr;
00452 char s_mo[3], s_da[3];
00453
00454
00455
00456
00457
00458
00459
00460
00461 if (SimTime.mo[0]<10) sprintf(s_mo,"0%d", SimTime.mo[0]);
00462 else sprintf(s_mo,"%d", SimTime.mo[0]);
00463 if (SimTime.da[0]<10) sprintf(s_da,"0%d", SimTime.da[0]);
00464 else sprintf(s_da,"%d", SimTime.da[0]);
00465 sprintf(ftype,"%d%s%s\0",SimTime.yr[0],s_mo,s_da);
00466
00467 if( HDF && type=='H' ) sprintf(ftype,".hdf");
00468
00469 gSize = gridSize*bSize + 1;
00470 if( gSize > gTempSize ) {
00471 if(gTempSize) free((char*)gTemp);
00472 gTemp = (UCHAR*) nalloc( gSize, "gTemp" );
00473 gTempSize = gSize;
00474 }
00475
00476 if( H_OPSYS == UNIX ){
00477 if(check_for((char*)filename, "/", 0, CASE, END ) >= 0 ) {
00478 if(filename[0] != '/' ) pathType = 1;
00479 else pathType = 2;
00480 }
00481 }
00482 else { if(check_for ((char*)filename, ":", 0, CASE, END ) >= 0 ) { if( filename[0] == ':' ) pathType = 1; else pathType = 2; }}
00483
00484 if(type=='H') {
00485 switch(pathType) {
00486 case 0:
00487 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Animation/%s%s\0",OutputPath,ProjName,filename,ftype);
00488 else sprintf(modelFileName,"%s%s:Output:Animation:%s%s\0",OutputPath,ProjName,filename,ftype);
00489 break;
00490 case 1:
00491 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Animation/%s%s\0",OutputPath,ProjName,filename,ftype);
00492 else sprintf(modelFileName,"%s%s:Output:Animation:%s%s\0",OutputPath,ProjName,filename,ftype);
00493 break;
00494 case 2:
00495 sprintf(modelFileName,"%s.hdf\0",filename);
00496 break;
00497 }
00498 }
00499 else if(type=='M') {
00500 switch(pathType) {
00501 case 0:
00502 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Maps/%s%s\0",OutputPath,ProjName,filename,ftype);
00503 else sprintf(modelFileName,"%s%s:Output:Maps:%s%s\0",OutputPath,ProjName,filename,ftype);
00504 break;
00505 case 1:
00506 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Maps/%s%s\0",OutputPath,ProjName,filename,ftype);
00507 else sprintf(modelFileName,"%s%s:Output:Maps:%s%s\0",OutputPath,ProjName,filename,ftype);
00508 break;
00509 case 2:
00510 sprintf(modelFileName,"%s%s\0",filename,ftype);
00511 break;
00512 }
00513 }
00514 else {
00515 bSize = 1;
00516 switch(pathType) {
00517 case 0:
00518 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Output/Animation%d/%s%s\0",OutputPath,ProjName,type,filename,ftype);
00519 else sprintf(modelFileName,"%s%s:Output:Animation%d:%s%s\0",OutputPath,ProjName,type,filename,ftype);
00520 break;
00521 case 1:
00522 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/DriverOutput/Animation%d/%s%s\0",OutputPath,ProjName,type,filename,ftype);
00523 else sprintf(modelFileName,"%s%s:Output:Animation%d:%s%s\0",OutputPath,ProjName,type,filename,ftype);
00524 break;
00525 case 2:
00526 sprintf(modelFileName,"%s%s\0",filename,ftype);
00527 break;
00528 }
00529 }
00530
00531 bSize = (bSize < 1) ? 1 : bSize;
00532 bSize = (bSize > 4) ? 4 : bSize;
00533 if( Mtype == 'c' ) bSize = 1;
00534 if(bSize == 1) { imax = 255; imin = 0; }
00535 else if(bSize == 2) { imax = 255*128; imin = -imax; }
00536 else if(bSize == 3) { imax = 256*256*128; imin = -imax; }
00537 mPtr = gTemp;
00538
00539 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",
00540 filename,scale_value,offset_value,index,modelFileName,Mtype,bSize,pathType,type); WriteMsg(msgStr,1); }
00541
00542 for (i=0; i<s0; i++) {
00543 for (j=0; j<s1; j++) {
00544 if( ON_MAP[T(i+1,j+1)] ) {
00545 switch(Mtype) {
00546 case 'l' :
00547 ftmp = ((( (double*) Map)[T(i+1,j+1)] - offset_value) / scale_value);
00548 itmp = (int) ftmp;
00549 itmp = (itmp > imax-1) ? imax-1 : itmp;
00550 itmp = (itmp < imin) ? imin : itmp;
00551 break;
00552 case 'f' :
00553 ftmp = ((( (float*) Map)[T(i+1,j+1)] - offset_value) / scale_value);
00554 itmp = (int) ftmp;
00555 itmp = (itmp > imax-1) ? imax-1 : itmp;
00556 itmp = (itmp < imin) ? imin : itmp;
00557 break;
00558 case 'd' : case 'i' :
00559 ftmp = ((( (int*) Map)[T(i+1,j+1)] - offset_value) / scale_value);
00560 itmp = (int) ftmp;
00561 itmp = (itmp > imax-1) ? imax-1 : itmp;
00562 itmp = (itmp < imin) ? imin : itmp;
00563 break;
00564 case 'c' :
00565 itmp = ((unsigned char*) Map)[T(i+1,j+1)];
00566 itmp = (itmp > imax-1) ? imax-1 : itmp;
00567 itmp = (itmp < imin) ? imin : itmp;
00568 break;
00569 }
00570 }
00571 else itmp = imax;
00572
00573 enc_Nb(&mPtr,itmp,bSize);
00574 }
00575 }
00576 writeMap(modelFileName,gTemp,bSize,type,index);
00577 if(debug) quick_look(Map, (char*)filename, Mtype, dbgPt.x, dbgPt.y, 2,'E');
00578 return(0);
00579 }
00580
00581
00591 int read_map_file(const char *filename, VOIDP Map,
00592 unsigned char Mtype, float scale_value,
00593 float offset_value)
00594 {
00595 int i, j, k0, size;
00596 unsigned temp;
00597 UCHAR *tmpPtr, Dsize;
00598
00599 gridSize = (s0+2)*(s1+2);
00600 gTempSize = gridSize*8;
00601
00602
00603 size = gridSize*4 +1;
00604 if( size > gTempSize ) {
00605 if(gTemp) free((char*)gTemp);
00606 gTemp = (UCHAR*)nalloc( size, "gTemp" );
00607 gTempSize = size;
00608 }
00609
00610 if(debug>2) {
00611 sprintf(msgStr,"Reading %s\n",filename);
00612 usrErr(msgStr);
00613 WriteMsg(msgStr,1);
00614 }
00615
00616 Dsize = readMap(filename,gTemp);
00617
00618 if(debug) {
00619 sprintf(msgStr,"name = %s, proc = %d, scale_value = %f, offset_value = %f, size = %x\n ",
00620 filename,procnum,scale_value,offset_value,Dsize);
00621 WriteMsg(msgStr,1);
00622 }
00623 for (i=1; i<=s0; i++) {
00624 for (j=1; j<=s1; j++) {
00625 k0 = Dsize*((i-1)*lcl_size[1] + (j-1));
00626 tmpPtr = gTemp+k0;
00627 switch(Dsize) {
00628 case 1: temp = gTemp[k0]; break;
00629 case 2: temp = gTemp[k0]*256 + gTemp[k0+1]; break;
00630 case 3: temp = gTemp[k0]*65536 + gTemp[k0+1]*256 + gTemp[k0+2]; break;
00631 case 4: temp = gTemp[k0]*16777216 + gTemp[k0+1]*65536 + gTemp[k0+2]*256 + gTemp[k0+3]; break;
00632 default: fprintf(stderr,"ERROR, illegal size: %x\n",Dsize); temp = 0;
00633 }
00634 switch (Mtype) {
00635 case 'f' :
00636 ((float*)Map)[T(i,j)] = scale_value*((float)temp)+offset_value;
00637 break;
00638 case 'd' : case 'i' :
00639 ((int*)Map)[T(i,j)] = (int)(scale_value * (float)temp + offset_value);
00640 break;
00641 case 'c' :
00642 ((unsigned char*)Map)[T(i,j)] = (int)(scale_value * (unsigned char)temp);
00643 break;
00644 }
00645 }
00646 }
00647 if(debug) quick_look(Map, (char*)filename, Mtype, dbgPt.x, dbgPt.y, 2,'E');
00648 link_edges(Map,Mtype);
00649 fflush(stderr);
00650 return Dsize;
00651 }
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00678 void enc_Nb(UCHAR** sptr,SLONG value,int bytes)
00679 {
00680 UCHAR tmp[4];
00681 switch(bytes) {
00682 case 1:
00683 *(*sptr)++ = value;
00684 break;
00685 case 2:
00686 *(*sptr)++ = (value >> 8);
00687 *(*sptr)++ = value;
00688 break;
00689 case 3:
00690 tmp[0] = value;
00691 tmp[1] = (value >>= 8);
00692 tmp[2] = (value >> 8);
00693 *(*sptr)++ = tmp[2];
00694 *(*sptr)++ = tmp[1];
00695 *(*sptr)++ = tmp[0];
00696 break;
00697 case 4:
00698 tmp[0] = value;
00699 tmp[1] = (value >>= 8);
00700 tmp[2] = (value >>= 8);
00701 tmp[3] = (value >> 8);
00702 *(*sptr)++ = tmp[3];
00703 *(*sptr)++ = tmp[2];
00704 *(*sptr)++ = tmp[1];
00705 *(*sptr)++ = tmp[0];
00706 break;
00707 }
00708 }
00709
00719 void quick_look(void* Map, char* name, unsigned char Mtype,
00720 int xc, int yc, int range, unsigned char format)
00721 {
00722 int ymin, ymax, xmin, xmax, x, y, N0, N1;
00723 char* namePtr;
00724
00725 if (!on_this_proc(xc,yc)) return;
00726 x = xc - lcl_start[0] + 1;
00727 y = yc - lcl_start[1] + 1;
00728 xmin = x - range;
00729 xmax = x + range + 1;
00730 ymin = y - range;
00731 ymax = y + range + 1;
00732 xmax = iclip( xmax, 0, s0+2);
00733 xmin = iclip( xmin, 0, s0+2);
00734 ymax = iclip( ymax, 0, s1+2);
00735 ymin = iclip( ymin, 0, s1+2);
00736 N0 = xmax-xmin;
00737 N1 = ymax-ymin;
00738 namePtr = name_from_path(name);
00739 if (debug >3) {
00740 sprintf(msgStr,"Window to Map %s, s=%d, proc0 = %d, proc1 = %d; Corners(gbl): (%d,%d) (%d,%d)\n",
00741 namePtr,range,recpnum[0],recpnum[1],xc-range,yc-range,xc+range+1,yc+range+1);
00742 writeWindow(Map,namePtr,msgStr,xmin,ymin,N0,N1,istep,Mtype,format);
00743 }
00744
00745 }
00746
00747
00760 void writeWindow(void* fValue, char* vName, char* desc,
00761 int x0, int y0, int N0, int N1, int index,
00762 UCHAR Mtype, UCHAR format)
00763 {
00764 int mSize, size;
00765 char ctemp[200];
00766
00767 switch(Mtype) {
00768 case 'l': size = sizeof(double); break;
00769 case 'f': case 'L': case 'E': size = sizeof(float); break;
00770 case 'd': case 'i': size = sizeof(int); break;
00771 case 'c': size = sizeof(char); break;
00772 }
00773 if( (mSize=size*N0*N1) > gTempSize ) {
00774 if(gTemp != NULL) free((char*)gTemp);
00775 gTemp = (UCHAR*)malloc(gTempSize=mSize);
00776 }
00777
00778 Copy(((UCHAR*)fValue)+T(x0,y0)*size, gTemp, N1*size, N0, (s1+2)*size, N1*size);
00779 sprintf(ctemp,"(%d)WIND: %s",index,vName);
00780 writeSeries(gTemp,ctemp,desc,N0,N1,Mtype,format);
00781 }
00782
00789 int iclip ( int x0, int min, int max )
00790 {
00791 int rv;
00792 rv = ( x0 > max) ? max : x0;
00793 rv = ( rv < min) ? min : rv;
00794 return rv;
00795 }
00796
00797
00806 void Copy(void *src, void *dst, int w, int n, int sw, int dw)
00807 {
00808 int i;
00809 for(i=0; i<n; i++)
00810 memcpy( ((UCHAR*)dst)+i*dw, ((UCHAR*)src)+i*sw, w );
00811 }
00812
00813
00814
00815
00817 void set_env_vars(void) {
00818
00819 FILE *infile;
00820 char filename[100], ch;
00821 static long start;
00822 int i, maxLen =200;
00823
00824 ModelPath = getenv("ModelPath");
00825 ProjName = getenv("ProjName");
00826 DriverPath = getenv("DriverPath");
00827 OS_TYPE = getenv("OSTYPE");
00828
00829 sprintf(msgStr,"OS type = %s ",OS_TYPE);
00830 usrErr(msgStr);
00831
00832 sprintf(msgStr,"Project Name = %s ",ProjName);
00833 usrErr(msgStr);
00834
00835 }
00836
00837
00846 void print_point(ViewParm *vp, void* Map, char Mtype,
00847 char* mName, int tIndex, int vpindex)
00848 {
00849 int x, y;
00850 Point3D *pt;
00851
00852 pt = vp->points + vpindex;
00853 x = pt->x - lcl_start[0];
00854 y = pt->y - lcl_start[1];
00855 if( (x >= 0) && (x < s0) && (y >= 0) && (y < s1) ) {
00856 if(tIndex==0) { pt->z = PListIndex; pSeries[PListIndex].Length=0;
00857 }
00858 else PListIndex = pt->z;
00859 if(tIndex==0) {
00860
00861 if(PListIndex >= MAX_PTSERIES ) fatal("Too many Point series.");
00862 else if( pSeries[PListIndex].data == NULL ) pSeries[PListIndex].data = (float*) malloc( (N_iter*sizeof(float)/vp->step)+2 );
00863 if( pSeries[PListIndex].data == NULL ) usrErr("out of Memory for Spatial Timeseries.");
00864 pSeries[PListIndex].Loc.x = pt->x;
00865 pSeries[PListIndex].Loc.y = pt->y;
00866 pSeries[PListIndex].outstep = vp->step;
00867 strcpy(pSeries[PListIndex].name,mName);
00868 if(debug >2) {
00869 sprintf(msgStr,"\nSetup Pointlist %d for %s(%d), step=%d, point=(%d,%d)\n", PListIndex, mName, vpindex, vp->step, x, y );
00870 WriteMsg(msgStr,1);
00871 }
00872 }
00873 switch(Mtype) {
00874 case 'f' : pSeries[PListIndex].data[ pSeries[PListIndex].Length++ ] = ((float*)Map)[T(x+1,y+1)]; break;
00875 case 'i' : pSeries[PListIndex].data[ pSeries[PListIndex].Length++ ] = ((int*)Map)[T(x+1,y+1)]; break;
00876 case 'c' : pSeries[PListIndex].data[ pSeries[PListIndex].Length++ ] = ((unsigned char*)Map)[T(x+1,y+1)];
00877 }
00878 if(debug >3) {
00879 sprintf(msgStr,"\nWriting Point %d for %s(%d), point=(%d,%d), value = %f, index = %d\n",
00880 PListIndex, mName, vpindex, x, y, ((float*)Map)[T(x+1,y+1)], pSeries[PListIndex].Length-1 );
00881 WriteMsg(msgStr,1);
00882 }
00883 PListIndex++;
00884 }
00885 }
00886
00892 void read_map_dims(const char* filename)
00893 {
00894 FILE *file;
00895
00896 if (debug>3) { sprintf(msgStr,"Getting map dims: %s",filename); usrErr(msgStr);}
00897 if(Lprocnum == 1) {
00898 sprintf(modelFileName,"%s/%s/Data/Map_head/%s",ModelPath,ProjName,filename);
00899 file = fopen(modelFileName,"r");
00900 if(file==NULL) {
00901 fprintf(stderr,"Unable to open map header file %s.\n",modelFileName);
00902 fflush(stderr);
00903 Exit(0);
00904 }
00905 scan_forward(file,"cols:");
00906 fscanf(file,"%d",&gbl_size[1]);
00907 sprintf(msgStr,"cols = %d\n",gbl_size[1]); WriteMsg(msgStr,1);
00908 scan_forward(file,"rows:");
00909 fscanf(file,"%d",&gbl_size[0]);
00910 sprintf(msgStr,"rows = %d\n",gbl_size[0]); WriteMsg(msgStr,1);
00911 fclose(file);
00912 }
00913 broadcastInt(gbl_size);
00914 broadcastInt(gbl_size+1);
00915 s0 = gbl_size[0];
00916 s1 = gbl_size[1];
00917 }
00918
00925 void init_pvar(VOIDP Map, UCHAR* mask, unsigned char Mtype,float iv)
00926 {
00927 int i0, i1;
00928
00929 switch(Mtype) {
00930 case 'b' :
00931 for(i0=0; i0<=numBasn; i0++) {
00932 ((double*)Map)[i0] = iv;
00933 }
00934 break;
00935 case 'l' :
00936 for(i0=0; i0<=s0+1; i0++)
00937 for(i1=0; i1<=s1+1; i1++) {
00938 if(mask==NULL) ((double*)Map)[T(i0,i1)] = iv;
00939 else if ( mask[T(i0,i1)] == 0 ) ((double*)Map)[T(i0,i1)] = 0;
00940 else ((double*)Map)[T(i0,i1)] = iv;
00941 }
00942 break;
00943 case 'f' :
00944 for(i0=0; i0<=s0+1; i0++)
00945 for(i1=0; i1<=s1+1; i1++) {
00946 if(mask==NULL) ((float*)Map)[T(i0,i1)] = iv;
00947 else if ( mask[T(i0,i1)] == 0 ) ((float*)Map)[T(i0,i1)] = 0;
00948 else ((float*)Map)[T(i0,i1)] = iv;
00949 }
00950 break;
00951 case 'i' : case 'd' :
00952 for(i0=0; i0<=s0+1; i0++)
00953 for(i1=0; i1<=s1+1; i1++) {
00954 if(mask==NULL) ((int*)Map)[T(i0,i1)] = (int)iv;
00955 else if ( mask[T(i0,i1)] == 0 ) ((int*)Map)[T(i0,i1)] = 0;
00956 else ((int*)Map)[T(i0,i1)] = (int)iv;
00957 }
00958 break;
00959 case 'c' :
00960 for(i0=0; i0<=s0+1; i0++)
00961 for(i1=0; i1<=s1+1; i1++) {
00962 if(mask==NULL) ((unsigned char*)Map)[T(i0,i1)] = (unsigned char) '\0' + (int) iv;
00963 else if ( mask[T(i0,i1)] == 0 ) ((unsigned char*)Map)[T(i0,i1)] = '\0';
00964 else ((unsigned char*)Map)[T(i0,i1)] = (unsigned char) '\0' + (int) iv;
00965 }
00966 break;
00967
00968 default :
00969 printf(" in default case\n");
00970 break;
00971 }
00972 }
00973
00982 void print_loc_ave(Point3D *vt, void* Map, char Mtype, char* mName, int tIndex)
00983 {
00984 int ix, iy, x, y, x0, y0, range, xmax, ymax, xmin, ymin, type = 11;
00985 float vout[4];
00986 double sum = 0, normCoef = 0;
00987
00988 if(Mtype != 'f') { fprintf(stderr,"Warning: Wrong type in Loc Ave: %s(%c)\n",mName,Mtype); return; }
00989 x0 = vt->x;
00990 y0 = vt->y;
00991 range = vt->z;
00992 x = x0 - lcl_start[0];
00993 y = y0 - lcl_start[1];
00994 xmin = x - range;
00995 xmax = x + range;
00996 ymin = y - range;
00997 ymax = y + range;
00998 xmax = (xmax > (s0)) ? (s0) : xmax;
00999 xmin = (xmin < 0 ) ? 0 : xmin;
01000 ymax = (ymax > (s1)) ? (s1) : ymax;
01001 ymin = (ymin < 0 ) ? 0 : ymin;
01002
01003 for(ix=xmin; ix<xmax; ix++) {
01004 for(iy=ymin; iy<ymax; iy++) {
01005 if( ON_MAP[T(ix+1,iy+1)] ) sum += ((float*)Map)[T(ix+1,iy+1)];
01006 normCoef += 1.0;
01007 }
01008 }
01009
01010 vout[0] = sum; vout[1] = normCoef;
01011
01012 switch (vt->type) {
01013 case 's' : Combine(vout,mName,1,kSUM,tIndex); break;
01014 case 'a' : Combine(vout,mName,2,kAVE,tIndex); break;
01015 case 'A' : Combine(vout,mName,2,kAVECUM,tIndex); break;
01016 case 'S' : Combine(vout,mName,1,kSUMCUM,tIndex); break;
01017 }
01018 }
01019
01020
01022 void PTS_SetFields (PTSeries* thisStruct, int ix, int iy,
01023 int index, int format, int col)
01024 {
01025 float ptstep = 1.0, *value;
01026 int n0;
01027
01028 thisStruct->fX = ix;
01029 thisStruct->fY = iy;
01030 if(thisStruct->fValue) free((char*)thisStruct->fValue);
01031
01032 thisStruct->fValue = readSeriesCol(modelFileName,format,index,&n0,&ptstep,col);
01033
01034 gNPtTs = thisStruct->fNpt = n0;
01035 gTS = thisStruct->fTS = ptstep;
01036 value = thisStruct->fValue;
01037 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); }
01038 }
01039
01041 void PTS_Free(PTSeries* thisStruct)
01042 {
01043 if(thisStruct->fValue) free((char*)thisStruct->fValue);
01044 }
01045
01046
01048 void PTS_CopyFields (PTSeries* thisStruct, PTSeries pV)
01049 {
01050 thisStruct->fX = pV.fX;
01051 thisStruct->fY = pV.fY;
01052 thisStruct->fValue = pV.fValue;
01053 }
01054
01056 PTSeriesList* PTSL_Init(int nSlots, int format)
01057 {
01058 int i;
01059 PTSeriesList* thisStruct;
01060 thisStruct = (PTSeriesList*) malloc(sizeof(PTSeriesList));
01061 thisStruct->fList = (PTSeries*) malloc(sizeof(PTSeries)*nSlots);
01062 for(i=0; i<nSlots; i++) thisStruct->fList[i].fValue = NULL;
01063 thisStruct->fNSlots = nSlots;
01064 thisStruct->fNSlotsUsed = 0;
01065 thisStruct->fNptTS = 0;
01066 thisStruct->fFormat = format;
01067 return(thisStruct);
01068 }
01069
01071 void PTSL_AddpTSeries(PTSeriesList* thisStruct, int x, int y, int index, int seriesNum, int col )
01072 {
01073 PTSeries* valListtemp;
01074 int i, newSlots;
01075
01076 if( seriesNum >= thisStruct->fNSlots ) {
01077 newSlots = Max(20,(int)(thisStruct->fNSlots * 0.2));
01078 thisStruct->fNSlots += newSlots;
01079 valListtemp = (PTSeries*) malloc(sizeof(PTSeries)*(thisStruct->fNSlots));
01080 for(i=0; i<thisStruct->fNSlotsUsed; i++) {
01081 valListtemp[i].fNpt = thisStruct->fNptTS;
01082 PTS_CopyFields(valListtemp+i,thisStruct->fList[i]);
01083 }
01084 for(i=thisStruct->fNSlotsUsed; i<thisStruct->fNSlots; i++) valListtemp[i].fValue = NULL;
01085 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01086 thisStruct->fList = valListtemp;
01087 }
01088
01089 PTS_SetFields( thisStruct->fList + seriesNum, x, y, index, thisStruct->fFormat, col );
01090 thisStruct->fNSlotsUsed = seriesNum;
01091 }
01092
01094 void PTSL_Free(PTSeriesList* thisStruct)
01095 {
01096 int i;
01097 if(thisStruct == NULL) return;
01098 for(i=0; i<thisStruct->fNSlotsUsed; i++) PTS_Free( thisStruct->fList +i );
01099 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01100 free((char*)thisStruct);
01101 thisStruct = NULL;
01102 }
01103
01105 float PTSL_GetInterpolatedValue0(PTSeriesList* thisStruct, int x, int y, int step)
01106 {
01107 int i, dx, dy, my_debug;
01108 PTSeries pV;
01109 float wpow;
01110 float weight, InterpValue=0.0, distance=0.0;
01111 int r;
01112
01113 wpow = GP_IDW_pow;
01114
01115
01116 for(i=0; i<thisStruct->fNSlotsUsed; i++) {
01117 pV = thisStruct->fList[i];
01118 dx = (pV.fX-(x+lcl_start[0]));
01119 dy = (pV.fY-(y+lcl_start[1]));
01120 r = (dx*dx + dy*dy);
01121 if( r == 0 ) return pV.fValue[step];
01122 weight = (wpow == 1.0) ? (1.0)/r : 1.0/(r*r);
01123
01124
01125
01126 InterpValue += pV.fValue[step]*weight;
01127 distance += weight;
01128 }
01129 if (distance > 0) InterpValue /= distance;
01130 else InterpValue = (thisStruct->fList[0]).fValue[step];
01131 return (float) InterpValue;
01132 }
01133
01135 void PTSL_ReadLists(PTSeriesList* thisStruct, const char *ptsFileName,
01136 int index, float* timeStep, int* nPtTS, int col )
01137 {
01138 char pathname[150], infilename[60], ss[201], ret = '\n';
01139 FILE* cfile;
01140 int ix, iy, tst, sCnt=0;
01141 if( Lprocnum == 1 ) {
01142 if(H_OPSYS==UNIX) sprintf(modelFileName,"%s/%s/Data/%s.pts",ModelPath,ProjName,ptsFileName);
01143 else sprintf(modelFileName,"%s%s:Data:%s.pts",ModelPath,ProjName,ptsFileName);
01144 cfile = fopen(modelFileName,"r");
01145 if(cfile==NULL) {fprintf(stdout,"\nERROR: Unable to open timeseries file %s\n",modelFileName); Exit(1); }
01146 else { sprintf(msgStr,"\nReading file %s\n",modelFileName); WriteMsg(msgStr,1); }
01147 if (debug > 2) {sprintf(msgStr,"Reading %s timeseries, column %d",modelFileName, col); usrErr(msgStr);}
01148
01149
01150
01151 fgets(ss,200,cfile);
01152 }
01153 while(1) {
01154 if( Lprocnum == 1 ) {
01155 tst = fscanf(cfile,"%d",&ix);
01156 if( tst > 0 ) {
01157 fscanf(cfile,"%d",&iy);
01158 tst = fscanf(cfile,"%s",infilename);
01159 sprintf(modelFileName,"%s/%s/Data/%s",ModelPath,ProjName,infilename);
01160 }
01161 }
01162 broadcastInt(&tst);
01163 if(tst<1) break;
01164 broadcastInt(&ix);
01165 broadcastInt(&iy);
01166 PTSL_AddpTSeries(thisStruct, ix, iy, index, sCnt, col);
01167 sCnt++;
01168 }
01169 if( Lprocnum == 1 ) fclose(cfile);
01170 *nPtTS = gNPtTs;
01171 *timeStep = gTS;
01172 }
01173
01175 void PTSL_CreatePointMap(PTSeriesList* pList,void* Map, unsigned char Mtype,
01176 int step, float scale)
01177 {
01178 int i0, i1;
01179
01180 switch(Mtype) {
01181 case 'f' :
01182 for(i0=0; i0<=s0+1; i0++)
01183 for(i1=0; i1<=s1+1; i1++)
01184 ((float*)Map)[T(i0,i1)] = (ON_MAP[T(i0,i1)]) ? PTSL_GetInterpolatedValue0(pList,i0,i1,step)*scale : 0.0 ;
01185 break;
01186 case 'i' : case 'd' :
01187 for(i0=0; i0<=s0+1; i0++)
01188 for(i1=0; i1<=s1+1; i1++)
01189 ((int*)Map)[T(i0,i1)] = (int) ( (ON_MAP[T(i0,i1)]) ? PTSL_GetInterpolatedValue0(pList,i0,i1,step)*scale : 0 );
01190 break;
01191 case 'c' :
01192 for(i0=0; i0<=s0+1; i0++)
01193 for(i1=0; i1<=s1+1; i1++)
01194 ((unsigned char*)Map)[T(i0,i1)] = (unsigned char) '\0' + (int) ( (ON_MAP[T(i0,i1)]) ? PTSL_GetInterpolatedValue0(pList,i0,i1,step)*scale : 0 );
01195 break;
01196 }
01197 }
01198
01200 float PTSL_GetInterpolatedValue1(PTSeriesList* thisStruct, int x, int y, int step)
01201 {
01202 int i, dx, dy;
01203 PTSeries pV;
01204 RPointList *pList;
01205 float wpow;
01206 double weight, InterpValue=0.0, distance=0.0;
01207 long r;
01208
01209 pList = RPL_Init( 20 );
01210 wpow = GP_IDW_pow;
01211 for(i=0; i<thisStruct->fNSlotsUsed; i++) {
01212 pV = thisStruct->fList[i];
01213 dx = (pV.fX-x);
01214 dy = (pV.fY-y);
01215 r = dx*dx + dy*dy;
01216 if( r == 0 ) return pV.fValue[step];
01217 RPL_AddrPoint(pList, dx, dy, r, pV.fValue[step]);
01218 weight = (wpow == 1.0) ? ((double)1.0)/r : pow(r,-wpow);
01219 InterpValue += pV.fValue[step]*weight;
01220 distance += weight;
01221 }
01222 RPL_Sort(pList);
01223 RPL_Free(pList);
01224 return (float) InterpValue;
01225 }
01226
01228 void RP_SetFields (RPoint* thisStruct, int ix, int iy, float r, float value)
01229 {
01230 thisStruct->fX = ix;
01231 thisStruct->fY = iy;
01232 thisStruct->fr = r;
01233 thisStruct->fValue = value;
01234 }
01235
01237 void RP_CopyFields (RPoint* thisStruct, RPoint* that)
01238 {
01239 thisStruct->fX = that->fX;
01240 thisStruct->fY = that->fY;
01241 thisStruct->fr = that->fr;
01242 thisStruct->fValue = that->fValue;
01243 }
01244
01246 void RP_SwapFields (RPoint* thisStruct, RPoint* that)
01247 {
01248 RPoint temp;
01249 RP_CopyFields (&temp, thisStruct);
01250 RP_CopyFields (thisStruct, that);
01251 RP_CopyFields (that, &temp);
01252 }
01253
01255 RPointList * RPL_Init(int nPoints)
01256 {
01257 RPointList *thisStruct;
01258 thisStruct = (RPointList*) malloc(sizeof(RPointList));
01259 thisStruct->fList = (RPoint*) malloc(sizeof(RPoint)*nPoints);
01260 thisStruct->fNPt = nPoints;
01261 thisStruct->fNPtUsed = 0;
01262 return thisStruct;
01263 }
01264
01266 void RPL_AddrPoint(RPointList *thisStruct, int x, int y, float r, float value)
01267 {
01268 RPoint* pointListtemp;
01269 int i, newPoints;
01270
01271 if( thisStruct->fNPtUsed >= thisStruct->fNPt ) {
01272 newPoints = 20;
01273 thisStruct->fNPt += newPoints;
01274 pointListtemp = (RPoint*) malloc(sizeof(RPoint)*(thisStruct->fNPt));
01275 for(i=0; i<thisStruct->fNPtUsed; i++) {
01276 RP_CopyFields(pointListtemp+i,thisStruct->fList+i);
01277 }
01278 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01279 thisStruct->fList = pointListtemp;
01280 }
01281
01282 RP_SetFields( thisStruct->fList + thisStruct->fNPtUsed, x, y, r, value);
01283 thisStruct->fNPtUsed++;
01284 }
01285
01287 void RPL_Sort(RPointList *thisStruct)
01288 {
01289 int i, go=1;
01290
01291 while(go) {
01292 go = 0;
01293 for(i=1; i<thisStruct->fNPtUsed; i++) if(thisStruct->fList[i-1].fr > thisStruct->fList[i].fr) {
01294 RP_SwapFields (thisStruct->fList+i-1, thisStruct->fList+i);
01295 go = 1;
01296 }
01297 }
01298 }
01299
01300
01302 void RPL_Free(RPointList* thisStruct)
01303 {
01304 if(thisStruct == NULL) return;
01305 if( thisStruct->fList != NULL ) free((char*)thisStruct->fList);
01306 free((char*)thisStruct);
01307 thisStruct = NULL;
01308 }
01309
01310
01322 void setFlag(ViewParm* vp, UINT mask, int value )
01323 {
01324 if(value) vp->flags |= mask; else vp->flags &= ~mask;
01325 }
01326
01328 int getFlag(ViewParm* vp, UINT mask)
01329 {
01330 if( vp->flags & mask ) return 1; else return 0;
01331 }
01332
01334 void setPrecision(ViewParm* vp, int value )
01335 {
01336 if(value/3) vp->flags |= PMASK2; else vp->flags &= ~PMASK2;
01337 if((value-1)%2) vp->flags |= PMASK1; else vp->flags &= ~PMASK1;
01338 }
01339
01341 int getPrecision(ViewParm* vp)
01342 {
01343 int rv = 1;
01344 if(vp->flags & PMASK2) rv +=2;
01345 if(vp->flags & PMASK1) rv +=1;
01346 return rv;
01347 }
01348
01360 int check_for(char *s0, const char *s1, int start, int cs, int rp)
01361 {
01362
01363
01364
01365
01366
01367 int i, j=0, k=-1, Len1 = strlen(s1), Len0 = strlen(s0);
01368 char t1, t2;
01369
01370 while(k<0) {
01371 k=0;
01372 for(i=start; i< Len0; ++i) {
01373 t1 = s0[i];
01374 t2 = s1[j];
01375 if(cs==NOCASE) { t1 = tolower(t1); t2 = tolower(t2); }
01376 if (t1 == t2) j++;
01377 else {j=0; k=0;}
01378 if(j==Len1) { k=1; break; }
01379 }
01380 }
01381 if(k<=0) return(-1);
01382 else if(rp==BEG) return(i-Len1+1);
01383 else return(i+1);
01384 }
01385
01386
01392 int isInteger (char *target_str) {
01393
01394 int i=-1,first_num=0;
01395 char ch;
01396
01397 while( (ch=target_str[++i]) != '\0' ) {
01398 if( isdigit(ch) ) first_num=1;
01399 if( (ch=='-' || ch=='+') && first_num ) return(0);
01400 if( !( isspace(ch) || isdigit(ch) || ch=='-' || ch=='+') ) return(0);
01401 }
01402 return(1);
01403 }
01404
01410 int isFloat (char *target_str) {
01411
01412 int i=-1,first_num=0;
01413 char ch;
01414
01415 while( (ch=target_str[++i]) != '\0' ) {
01416 if( isdigit(ch) ) first_num=1;
01417 if( (ch=='-' || ch=='+') && first_num ) return(0);
01418 if( !( isspace(ch) || isdigit(ch) || ch=='-' || ch=='+' || ch=='.' || toupper(ch) == 'E') ) return(0);
01419 }
01420 return(1);
01421 }
01422
01423
01436 void WriteMsg(const char* msg, int wh) {
01437 wh = 1;
01438 if(Driver_outfile) fprintf(Driver_outfile,"%s\n",msg);
01439 else fprintf(stdout,"%s\n",msg);
01440 fflush(stdout);
01441 }
01442
01447 void usrErr0(const char* dString)
01448 {
01449 fprintf(stderr,"%s", dString);
01450 fflush(stderr);
01451 }
01452
01457 void usrErr(const char* dString)
01458 {
01459 fprintf(stderr,"%s\n", dString);
01460 fflush(stderr);
01461 }
01462
01463
01465 void Exit(int code) { exit(code); }
01466
01468 void fatal(const char *msg)
01469 {
01470 printf("%s",msg);
01471 fflush(stdout);
01472 Exit(9);
01473 }
01474
01475
01476
01477
01478
01500 double julday(int mon, int day, int year, int h, int mi, double se)
01501 {
01502 long m = mon, d = day, y = year;
01503 long c, ya, j;
01504 double seconds = h * 3600.0 + mi * 60 + se;
01505
01506 if (m > 2)
01507 m -= 3;
01508 else {
01509 m += 9;
01510 --y;
01511 }
01512 c = y / 100L;
01513 ya = y - (100L * c);
01514 j = (146097L * c) / 4L + (1461L * ya) / 4L + (153L * m + 2L) / 5L + d + 1721119L;
01515 if (seconds < 12 * 3600.0) {
01516 j--;
01517 seconds += 12.0 * 3600.0;
01518 }
01519 else {
01520 seconds = seconds - 12.0 * 3600.0;
01521 }
01522 return (j + (seconds / 3600.0) / 24.0);
01523 }
01524
01547 void calcdate(double jd, int *m, int *d, int *y, int *h, int *mi, double *sec)
01548 {
01549 static int ret[4];
01550
01551 long j = (long)jd;
01552 double tmp, frac = jd - j;
01553
01554 if (frac >= 0.5) {
01555 frac = frac - 0.5;
01556 j++;
01557 }
01558 else {
01559 frac = frac + 0.5;
01560 }
01561
01562 ret[3] = (j + 1L) % 7L;
01563 j -= 1721119L;
01564 *y = (4L * j - 1L) / 146097L;
01565 j = 4L * j - 1L - 146097L * *y;
01566 *d = j / 4L;
01567 j = (4L * *d + 3L) / 1461L;
01568 *d = 4L * *d + 3L - 1461L * j;
01569 *d = (*d + 4L) / 4L;
01570 *m = (5L * *d - 3L) / 153L;
01571 *d = 5L * *d - 3 - 153L * *m;
01572 *d = (*d + 5L) / 5L;
01573 *y = 100L * *y + j;
01574 if (*m < 10)
01575 *m += 3;
01576 else {
01577 *m -= 9;
01578 *y += 1;
01579 }
01580 tmp = 3600.0 * (frac * 24.0);
01581 *h = (int) (tmp / 3600.0);
01582 tmp = tmp - *h * 3600.0;
01583
01584 *mi = (int) (tmp / 60.0);
01585 *sec = tmp - *mi * 60.0;
01586 }
01587
01592 float FMOD(float x, float y) {
01593 return (float)fmod((double)x, (double)y);
01594 }
01595
01600 void getFloat(FILE* inFile, const char* lString, float* fValPtr)
01601 {
01602 int test, fSize;
01603 UCHAR iEx=0;
01604 if(lString) scan_forward(inFile,lString);
01605 test = fscanf(inFile,"%f",fValPtr);
01606 if(test != 1) {
01607 fprintf(stderr,"Read Error (%d):",test);
01608 if(lString) fprintf(stderr,"%s\n",lString);
01609 *fValPtr = 0.0;
01610 iEx = 1;
01611 }
01612 if (iEx) exit(0);
01613 }
01614
01619 void getChar(FILE* inFile, const char* lString, char* cValPtr)
01620 {
01621 int test;
01622 UCHAR iEx=0;
01623 if(lString) scan_forward(inFile,lString);
01624 test = fscanf(inFile,"%c",cValPtr);
01625 if(test != 1) {
01626 fprintf(stderr,"Read Error (%d):",test);
01627 if(lString) fprintf(stderr,"%s\n",lString);
01628 iEx = 1;
01629 *cValPtr = '\0';
01630 }
01631 if (iEx) exit(0);
01632 }
01633
01634
01639 void getString(FILE* inFile, const char *lString, char *inString)
01640 {
01641 int test;
01642 UCHAR iEx=0;
01643 if(lString) scan_forward(inFile,lString);
01644 test = fscanf(inFile,"%s",inString);
01645 if(test != 1) {
01646 fprintf(stderr,"Read Error (%d):",test);
01647 if(lString) fprintf(stderr,"%s\n",lString);
01648 iEx = 1;
01649 }
01650 if (iEx) exit(0);
01651 }
01652
01657 void getInt(FILE* inFile, const char* lString, int* iValPtr)
01658 {
01659 int test;
01660 UCHAR iEx=0;
01661 if(lString)
01662 scan_forward(inFile,lString);
01663 test = fscanf(inFile,"%d",iValPtr);
01664 if(test != 1) {
01665 fprintf(stderr,"Read Error (%d):",test);
01666 if(lString) fprintf(stderr,"%s\n",lString);
01667 *iValPtr = 0;
01668 iEx = 1;
01669 }
01670 if (iEx) exit(0);
01671 }
01672
01677 int find_char(FILE *infile, char tchar)
01678 {
01679 char test = '1', cchar = '#';
01680 int in_c=0;
01681
01682 while( test != EOF ) {
01683 test = fgetc(infile);
01684 if(test == tchar && in_c==0 ) { return 1; }
01685 if(test == cchar ) in_c=1;
01686 if(test == '\n' ) in_c=0;
01687 }
01688 return(0);
01689 }
01690
01694 int Round(float x)
01695 {
01696 int i = (int)x;
01697 return (i);
01698 }
01699
01705 char *Scip(char *s, char SYM )
01706 {
01707 if(*s == SYM ) return ++s;
01708 while(*s != SYM && *s != EOS ) s++;
01709 if(*s != EOS) return ++s;
01710 else return s;
01711 }
01712
01713
01717 int skip_white(FILE *infile)
01718 {
01719 int ch;
01720
01721 while( isspace(ch=fgetc(infile)) ) {;}
01722 if(ch==EOF) return 0;
01723 ungetc(ch,infile);
01724 return 1;
01725 }
01726
01731 int scan_forward ( FILE *infile, const char *tstring)
01732 {
01733 int sLen, i, cnt=0;
01734 char Input_string[100], test;
01735
01736 sLen = strlen(tstring);
01737 while( ( test = fgetc(infile) ) != EOF ) {
01738 for(i=0; i<(sLen-1); i++)
01739 Input_string[i] = Input_string[i+1];
01740 Input_string[sLen-1] = test;
01741 Input_string[sLen] = '\0';
01742 if(++cnt >= sLen) {
01743 test = strcmp(Input_string,tstring);
01744 if( test == 0 ) return 1;
01745 }
01746 }
01747 return(-1);
01748 }
01749
01750
01754 char* name_from_path(char* name)
01755 {
01756 char* namePtr; int i, slen;
01757 char dirCh;
01758
01759 namePtr = name;
01760 dirCh = '/';
01761 slen = strlen(name);
01762
01763 for(i=0; i<slen; i++) {
01764 if( name[slen-i-1] == dirCh ) { namePtr = name+slen-i; break; }
01765 }
01766 return namePtr;
01767 }
01768
01769
01774 VOIDP nalloc(unsigned mem_size, const char var_name[])
01775 {
01776 VOIDP rp;
01777
01778
01779 if(mem_size == 0) return(NULL);
01780 rp = (VOIDP)malloc( mem_size );
01781 total_memory += mem_size;
01782 fasync(stderr);
01783 if( rp == NULL ) {
01784 fprintf(stderr,"Sorry, out of memory(%d): %s\n",mem_size,var_name);
01785 Exit(0);
01786 }
01787 fmulti(stderr);
01788 return(rp);
01789 }
01790
01792 float Normal(float mean,float sd)
01793 {
01794 int table_loc ;
01795 double rand_num ;
01796 float sign ;
01797 float interval = .1 ;
01798 int n_interval = 40 ;
01799 float high, low ;
01800 float return_val = 0.0;
01801
01802 sign = SMDRAND(0.0, 1.0) ;
01803 sign = (sign < .5 ? -1 : 1) ;
01804 rand_num = SMDRAND(.5, 1.0);
01805 low = gRTable[0] ;
01806 for (table_loc=1; table_loc<n_interval; table_loc++)
01807 {
01808 high = gRTable[table_loc] ;
01809 if (high > rand_num + .5)
01810 {
01811 return_val = mean + sd * (sign * interval * (table_loc - 1 +
01812 (rand_num+.5-low)/(high-low))) ;
01813 return(return_val) ;
01814 }
01815 low = high ;
01816 }
01817 return(return_val) ;
01818 }
01819
01821 int Poisson(float mu)
01822 {
01823 int ix;
01824 float f0, r, Lp = 1, xf = 1, p=0;
01825
01826 p = f0 = Exp(-mu);
01827 r = SMDRAND(0.0,1.0);
01828
01829 for( ix=0; ix<500; ix++) {
01830 if( r < p ) return ix;
01831 Lp *= mu;
01832 if(ix>0) xf *= ix;
01833 p += (Lp*f0)/xf;
01834 }
01835 return ix;
01836 }
01837
01838
01839
01841 void link_edges(VOIDP Map, unsigned char Mtype)
01842 {
01843 switch(Mtype) {
01844 case 'f' :
01845 exchange_borders((byte*)Map,sizeof(float)); break;
01846 case 'd' : case 'i' :
01847 exchange_borders((byte*)Map,sizeof(int)); break;
01848 case 'c' :
01849 exchange_borders((byte*)Map,sizeof(UCHAR)); break;
01850 }
01851 }
01852
01854 void setup_platform() {
01855
01856 exparam(&env);
01857 procnum = env.procnum;
01858 exgridsplit(env.nprocs, 2, nprocs);
01859 if( exgridinit(2, nprocs) < 0) { usrErr("Failed to Setup Grid"); exit(0); }
01860 exgridcoord(procnum, recpnum);
01861 }
01862
01863
01864
01865