void create_hfile(string hfile, string dfile, grid g, float val, bool var=false) {
        
        if (hfile.size()>120) {
            RVLException e;
            e<<"Error: create_hfile\n";
            e<<"filename "<<hfile<<" longer than 120 chars\n";
            throw e;
        }
        
        char * fname=(char *)malloc(128*sizeof(char));
        FILE * fp = NULL;
        
        strcpy(fname,hfile.c_str());
        
        fp = iwave_fopen(&fname,"w",NULL,stderr);
        if (!fp) {
            RVLException e;
            e<<"Error: create_hfile\n";
            e<<"file testgrid.rsf not opened\n";
            throw e;
        }
        
        fprint_grid(fp,g);
        
        fprintf(fp,"data_format=native_float\n");
        
        fprintf(fp,"data_type = csq\n");
        fprintf(fp,"in=%s\n",dfile.c_str());
        
        fflush(fp);
        
        iwave_fclose(fp);
        
        strcpy(fname,dfile.c_str());
        float * buf = (float *)malloc(get_datasize_grid(g)*sizeof(float));
        if (var) {
            for (int i=0;i<get_datasize_grid(g);i++) buf[i]=val/((float)(i+1));
        }
        else {
            for (int i=0;i<get_datasize_grid(g);i++) buf[i]=val;
        }
        //      FILE *
        fp = iwave_fopen(&fname,"w",NULL,stderr);
        if (!fp) {
            RVLException e;
            e<<"Error: create_hfile\n";
            e<<"file "<<fname<<" not opened\n";
            throw e;
        }
        for (int i=0;i<get_panelnum_grid(g);i++) 
            fwrite(buf,sizeof(float),get_datasize_grid(g),fp);
        
        free(buf);
        fflush(fp);
        
        iwave_fclose(fp);    
        free(fname);
    }

  void create_gauss(grid g, float lam) {
    
    char * fname=(char *)malloc(128*sizeof(char));
    FILE * fp = NULL;

    strcpy(fname,"gauss.rsf");
      
    fp = iwave_fopen(&fname,"w",NULL,stderr);
    if (!fp) {
      RVLException e;
      e<<"Error: create_hfile\n";
      e<<"file gauss.rsf not opened\n";
      throw e;
    }
    
    fprint_grid(fp,g);
    
    fprintf(fp,"data_format=native_float\n");
    
    fprintf(fp,"data_type = csq\n");
    fprintf(fp,"in=%s\n","gauss.rsf@");
    
    fflush(fp);
    
    iwave_fclose(fp);
    
    strcpy(fname,"gauss.rsf@");
    
    // compute center
    RPNT ctr;
    float radsq=0.0f;
    
    for (int i=0;i<g.dim;i++) {
      ctr[i] = g.axes[i].o + 0.5 * g.axes[i].d * g.axes[i].n;
      radsq += (0.5 * g.axes[i].d * g.axes[i].n)*(0.5 * g.axes[i].d * g.axes[i].n);
    }
    float x;
    float y; 
    float z;
    float rsq;
    float * buf = (float *)malloc(get_datasize_grid(g)*sizeof(float));
    if (g.dim==2) {
      //	cerr<<"g.axes[0].n="<<g.axes[0].n<<" g.axes[1].n="<<g.axes[1].n<<endl;
      //	cerr<<"g.axes[0].o="<<g.axes[0].o<<" g.axes[1].o="<<g.axes[1].o<<endl;
      //	cerr<<"g.axes[0].d="<<g.axes[0].d<<" g.axes[1].d="<<g.axes[1].d<<endl;
      //	cerr<<"radsq="<<radsq<<endl;
      for (int j=0;j<g.axes[1].n;j++) {
	for (int i=0;i<g.axes[0].n;i++) {
	  z=g.axes[0].o + i*g.axes[0].d;
	  x=g.axes[1].o + j*g.axes[1].d;
	  rsq=(z-ctr[0])*(z-ctr[0])+(x-ctr[1])*(x-ctr[1]);
	  //	    cerr<<"rsq/radsq="<<rsq/radsq<<endl;
	  buf[i+j*g.axes[0].n]=exp(-lam*lam*rsq/radsq);
	}
      }
    }
    else if (g.dim==3) {
      for (int k=0;k<g.axes[2].n;k++) {
	for (int j=0;j<g.axes[1].n;j++) {
	  for (int i=0;i<g.axes[0].n;i++) {
	    z=g.axes[0].o + i*g.axes[0].d;
	    x=g.axes[1].o + j*g.axes[1].d;
	    y=g.axes[2].o + k*g.axes[2].d;
	    buf[i+j*g.axes[0].n]=exp(-lam*lam*((z-ctr[0])*(z-ctr[0])+(x-ctr[1])*(x-ctr[1])+(y-ctr[2])*(y-ctr[2])));
	  }
	}
      }  
    }
    else {
      RVLException e;
      e<<"Error: create_hfile\n";
      e<<"  wake up and smell the roses\n";
      throw e;
    }
    //      FILE * 
    fp = iwave_fopen(&fname,"w",NULL,stderr);
    if (!fp) {
      RVLException e;
      e<<"Error: create_hfile\n";
      e<<"file "<<fname<<" not opened\n";
      throw e;
    }
    for (int i=0;i<get_panelnum_grid(g);i++) 
      fwrite(buf,sizeof(float),get_datasize_grid(g),fp);
    
    free(buf);
    fflush(fp);
    
    iwave_fclose(fp);
    free(fname);
  }