Everglades Landscape Model (ELM) Home Page
Main Page | Data Structures | Directories | File List | Data Fields | Globals

Driver_Utilities.c

Go to the documentation of this file.
00001 
00020 /* General notes on revisions to this source file. 
00021        Nov/Dec 2004 v2.3.2: documentation upgrade 
00022                 - Unused functions removed 
00023                 - Re-organized, clarified scopes, Doxygen tags added 
00024                 - Functionality same as v2.1/v2.2 
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 /* need to change the Outlist_size (#defined in header) if adding outlist variables to the data struct in gen_output of UnitMod.c */
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);/*fix*/
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               /* TODO: develop flexible output of calendar-based and julian-based outsteps (jan 11, 2005) */
00160               /* kluge here in checking use of Model.outList outstep relative to Driver.parm avg_Intvl */
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               /* TODO: develop flexible output of calendar-based and julian-based outsteps (jan 11, 2005) */
00221               /* kluge here in preventing use of point time series output using calendar-based outstep  */
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       } /* end of switch */
00282     } /* end of else */
00283   } /* end of while */
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); /* Iterative and cumulative, maximum and minimum values for all window cells at outstep iterations */
00373 }
00374 
00376 void setup_grid() {
00377   
00378   exgridsize(procnum,gbl_size,lcl_size,lcl_start); /* effectively unused in serial (non-parallel) implementation */
00379 
00380   /*  s0 and s1 remain unchanged in serial (non-parallel) implementation */
00381   s0 = lcl_size[0];
00382   s1 = lcl_size[1];
00383   
00384   gridSize = (s0+2)*(s1+2); /* the gridSize is increased by 2 in each dimension for buffer strips around processor domain(s) (functional mainly to parallel implementation) */
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[7] */ftype[10], gSize;
00449   float ftmp; 
00450   SLONG itmp, imax, imin;
00451   UCHAR *mPtr;
00452   char s_mo[3], s_da[3];
00453   
00454 /* below provides old (v1.0) increasing integer, non-date, appendage to filenames */
00455 /*   ftype[0] = '0' + (index / 100); */
00456 /*   ftype[1] = '0' + ((index%100) / 10); */
00457 /*   ftype[2] = '0' + (index % 10); */
00458  /*  ftype[8] = '\0'; */
00459 
00460 /*  append the date (yyyymmdd) to the root filename*/
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') {                         /* HDF format */
00485     switch(pathType) {
00486     case 0:     /* No Path */
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:     /* Relative Path */
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:     /* Full Path */
00495       sprintf(modelFileName,"%s.hdf\0",filename);
00496       break;
00497     }
00498   }
00499   else if(type=='M') {  
00500     switch(pathType) {                   /* Map II (float, binary) Format */
00501     case 0:                                                             /* No Path */
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:     /* Relative Path */
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:     /* Full Path */
00510       sprintf(modelFileName,"%s%s\0",filename,ftype);
00511       break;
00512     }
00513   }
00514   else {
00515         bSize = 1;
00516     switch(pathType) {                   /* Generic Binary Format */
00517     case 0:                                                             /* No Path */
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:     /* Relative Path */
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:     /* Full Path */
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); /* v2.3.2 added double (casting to float ftmp) */
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 /*      Routine:        enc_Nb                                  */
00658 /*                                                                      */
00659 /*      Place a 1,2,3, or 4 (N) byte signed value into the specified buffer     */
00660 /*      and advance the buffer pointer past the placed object.          */
00661 /*                                                                      */
00662 /*      Parameters:                                                     */
00663 /*              "sptr"  A pointer to pointer to UCHAR (the buffer).     */
00664 /*              "value" The value to place in the buffer.               */
00665 /*                                                                      */
00666 /*      Returns:        None.                                           */
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; /* v2.3.2 added */
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 ); /* copy w chars from source to dest, return destination */
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"); /* OSTYPE not used in code, only here for informational output */
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               /* TODO: develop flexible output of calendar-based and julian-based outsteps (jan 11, 2005) */
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);   /* does nothing in serial (non-parallel) implementation */
00914   broadcastInt(gbl_size+1); /* does nothing in serial (non-parallel) implementation */
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' :    /* added double for (non-map) basin (b) array budget calcs */
00931     for(i0=0; i0<=numBasn; i0++) {
00932         ((double*)Map)[i0] = iv;
00933       }
00934     break;
00935   case 'l' :    /* added double (l == letter "ell" ) for map arrays */
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; /* Sum for all window cells at outStep iterations */
01014      case 'a' :  Combine(vout,mName,2,kAVE,tIndex); break; /* Average for all window cells at outStep iterations */
01015      case 'A' :  Combine(vout,mName,2,kAVECUM,tIndex); break; /* Cumulative average for all window cells over all outStep iterations */
01016      case 'S' :  Combine(vout,mName,1,kSUMCUM,tIndex); break; /* Cumulative sum for all window cells over all outStep iterations */
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; /* was double in v1.0 */
01111   int r;
01112   
01113   wpow = GP_IDW_pow; /* GP_IDW_pow is "adjustable" parameter in global parm input file */
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         /* the pow() and double was incredibly slow (signif slowed simualtion), so removed it
01124            and gave the choice of only inverse distance or inverse distance**2 weighting */
01125     /* weight = (wpow == 1.0) ? ( (double)1.0)/r : pow(r,-wpow) ; */
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); /* skip header line */
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 /*was this value (v0.5beta) -0.012 */; 
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; /* GP_IDW_pow is "adjustable" parameter in global parm input file */
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   /*  Check for occurrences of string s1 in string s0   */
01363   /*  after position start. Return -1 if not found.             */
01364   /* if cs = CASE -> case sensitive, cs = NOCASE, not case sens.        */
01365   /* if rp = BEG -> return position of beginning of s1,         */
01366   /*                    otherwise return position of (next char following the) end of s1                */
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 /* The two julday & calcdate functions for going to/from julianDay/CalendarDate
01476    were (unmodified) from SFWMD HSM's /vol/hsm/src/libs/xmgr_julday/ directory.
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 

Generated on Thu Jul 6 11:17:44 2006 for ELM source code by  doxygen 1.3.9.1