//.........这里部分代码省略.........
	out_fd = G_open_raster_new(seedmap, 1);
	if (out_fd < 0)
	    G_fatal_error(_("Unable to create raster map <%s>"), seedmap);
    }

    /* More pases are renudant. Real pases count is controled by altered cell count. */
    pases = (int)(rows * cols) / 2;

    G_debug(1,
	    "Starting lake filling at level of %8.4f in %d passes. Percent done:",
	    water_level, pases);

    lastcount = 0;

    for (pass = 0; pass < pases; pass++) {
	G_debug(3, "Pass: %d", pass);
	curcount = 0;
	/* Move from left upper corner to right lower corner. */
	for (row = 0; row < rows; row++) {
	    for (col = 0; col < cols; col++) {
		/* Loading water data into window. */
		load_window_values(out_water, water_window, rows, cols, row,
				   col);

		/* Cheking presence of water. */
		if (is_near_water(water_window) == 1) {
		    if (in_terran[row][col] < water_level) {
			out_water[row][col] =
			    water_level - in_terran[row][col];
			curcount++;
		    }
		    else {
			out_water[row][col] = 0;	/* Cell is higher than water level -> NULL. */
		    }
		}
	    }
	}
	if (curcount == lastcount)
	    break;		/* We done. */
	lastcount = curcount;
	curcount = 0;
	/* Move backwards - from lower right corner to upper left corner. */
	for (row = rows - 1; row >= 0; row--) {
	    for (col = cols - 1; col >= 0; col--) {
		load_window_values(out_water, water_window, rows, cols, row,
				   col);

		if (is_near_water(water_window) == 1) {
		    if (in_terran[row][col] < water_level) {
			out_water[row][col] =
			    water_level - in_terran[row][col];
			curcount++;
		    }
		    else {
			out_water[row][col] = 0;
		    }
		}
	    }
	}
	G_percent(pass + 1, pases, 10);
	if (curcount == lastcount)
	    break;		/* We done. */
	lastcount = curcount;
    }				/*pases */

    G_percent(pases, pases, 10);	/* Show 100%. */

    save_map(out_water, out_fd, rows, cols, negative_flag->answer, &min_depth,
	     &max_depth, &area, &volume);

    G_message(_("Lake depth from %f to %f"), min_depth, max_depth);
    G_message(_("Lake area %f square meters"), area);
    G_message(_("Lake volume %f cubic meters"), volume);
    G_warning(_("Volume is correct only if lake depth (terrain raster map) is in meters"));

    /* Close all files. Lake map gets written only now. */
    G_close_cell(in_terran_fd);
    G_close_cell(out_fd);

    /* Add blue color gradient from light bank to dark depth */
    G_init_colors(&colr);
    if (negative_flag->answer == 1) {
	G_add_f_raster_color_rule(&max_depth, 0, 240, 255,
				  &min_depth, 0, 50, 170, &colr);
    }
    else {
	G_add_f_raster_color_rule(&min_depth, 0, 240, 255,
				  &max_depth, 0, 50, 170, &colr);
    }

    if (G_write_colors(lakemap, G_mapset(), &colr) != 1)
	G_fatal_error(_("Unable to read color file of raster map <%s>"),
		      lakemap);

    G_short_history(lakemap, "raster", &history);
    G_command_history(&history);
    G_write_history(lakemap, &history);

    return EXIT_SUCCESS;
}

void wtrshed(int fm, int fd, int nl, int ns, int mxbuf)
{
    int pass, repeat, flag, i, j, half, bufsz;
    int sline, nline, rdline;

    struct whereandwhat hold;
    struct whereandwhat *dir;
    struct whereandwhat *bas;

    dir = G_malloc(mxbuf * sizeof(struct whereandwhat));
    bas = G_malloc(mxbuf * sizeof(struct whereandwhat));

    bufsz = ns * sizeof(CELL);

    /* adjust maxbuf to an even number */
    half = mxbuf / 2;
    mxbuf = 2 * half;

    /* allocate buffers for drainage directions and basin areas */
    for (i = 0; i < mxbuf; i += 1)
        bas[i].p = (CELL *) G_calloc(ns, sizeof(CELL));
    for (i = 0; i < mxbuf; i += 1)
        dir[i].p = (CELL *) G_calloc(ns, sizeof(CELL));

    pass = 0;

    /* complete a downward pass */
    do {
        G_message(_("wtrshed pass %d"), ++pass);
        repeat = 0;

        /* fill the buffer */
        nline = mxbuf;
        sline = 0;
        rdline = 1;
        for (i = 0; i < mxbuf; i++) {
            bas[i].offset = dir[i].offset = (off_t) rdline *bufsz;

            lseek(fm, bas[i].offset, SEEK_SET);
            read(fm, bas[i].p, bufsz);

            lseek(fd, dir[i].offset, SEEK_SET);
            read(fd, dir[i].p, bufsz);

            rdline++;
        }

        /* repeat for all subsequent rows except the first and last */
        for (i = 1; i < nl - 1; i += 1) {
            /* analyse one line */
            for (j = 1; j < ns - 1; j += 1) {
                flag = bas[sline].p[j];
                if (flag > 0)
                    if (recurse_cell(flag, sline, j, nline, ns, bas, dir) > 0)
                        repeat = 1;
            }

            /* write one line */
            lseek(fm, bas[sline].offset, SEEK_SET);
            write(fm, bas[sline].p, bufsz);

            /* If the bottom end of the buffers reach the bottom of the file,
             * rotate the buffers and read new lines */
            if (rdline < nl - 1) {
                hold = bas[0];
                for (j = 1; j < mxbuf; j += 1)
                    bas[j - 1] = bas[j];
                bas[mxbuf - 1] = hold;

                hold = dir[0];
                for (j = 1; j < mxbuf; j += 1)
                    dir[j - 1] = dir[j];
                dir[mxbuf - 1] = hold;

                bas[mxbuf - 1].offset = dir[mxbuf - 1].offset =
                                            (off_t) rdline *bufsz;

                lseek(fm, bas[mxbuf - 1].offset, SEEK_SET);
                read(fm, bas[mxbuf - 1].p, bufsz);

                lseek(fd, dir[mxbuf - 1].offset, SEEK_SET);
                read(fd, dir[mxbuf - 1].p, bufsz);

                rdline++;
            }
            /* After the buffer reaches the bottom of the file, stop reading file,
             * just advance the pointers */
            else {
                nline -= 1;
                sline += 1;
            }

        }

        /* fill the buffer */
        nline = mxbuf;
        rdline = nl - 2;
        for (i = mxbuf - 1; i >= 0; i -= 1) {
            bas[i].offset = dir[i].offset = (off_t) rdline *bufsz;

//.........这里部分代码省略.........

//.........这里部分代码省略.........
	count = G_calloc(ncats, sizeof(DCELL));
	sum = G_calloc(ncats, sizeof(DCELL));
	break;
    case VARIANCE1:
    case STDDEV1:
	count = G_calloc(ncats, sizeof(DCELL));
	sum = G_calloc(ncats, sizeof(DCELL));
	sum2 = G_calloc(ncats, sizeof(DCELL));
	break;
    case SKEWNESS1:
	count = G_calloc(ncats, sizeof(DCELL));
	sum = G_calloc(ncats, sizeof(DCELL));
	sum2 = G_calloc(ncats, sizeof(DCELL));
	sum3 = G_calloc(ncats, sizeof(DCELL));
	break;
    case KURTOSIS1:
	count = G_calloc(ncats, sizeof(DCELL));
	sum = G_calloc(ncats, sizeof(DCELL));
	sum2 = G_calloc(ncats, sizeof(DCELL));
	sum4 = G_calloc(ncats, sizeof(DCELL));
	break;
    }

    if (min)
	for (i = 0; i < ncats; i++)
	    min[i] = 1e300;
    if (max)
	for (i = 0; i < ncats; i++)
	    max[i] = -1e300;

    base_buf = Rast_allocate_c_buf();
    cover_buf = Rast_allocate_d_buf();

    G_message(_("First pass"));

    for (row = 0; row < rows; row++) {
	Rast_get_c_row(base_fd, base_buf, row);
	Rast_get_d_row(cover_fd, cover_buf, row);

	for (col = 0; col < cols; col++) {
	    int n;
	    DCELL v;

	    if (Rast_is_c_null_value(&base_buf[col]))
		continue;
	    if (Rast_is_d_null_value(&cover_buf[col]))
		continue;

	    n = base_buf[col] - mincat;

	    if (n < 0 || n >= ncats)
		continue;

	    v = cover_buf[col];
	    if (usecats)
		sscanf(Rast_get_c_cat((CELL *) &v, &cats), "%lf", &v);

	    if (count)
		count[n]++;
	    if (sum)
		sum[n] += v;
	    if (sum2)
		sum2[n] += v * v;
	    if (sum3)
		sum3[n] += v * v * v;
	    if (sum4)

int execute_random(struct rr_state *theState)
{
    long nt;
    long nc;
    struct Cell_head window;
    int nrows, ncols, row, col;
    int infd, cinfd, outfd;
    struct Map_info Out;
    struct field_info *fi;
    dbTable *table;
    dbColumn *column;
    dbString sql;
    dbDriver *driver;
    struct line_pnts *Points;
    struct line_cats *Cats;
    int cat;
    RASTER_MAP_TYPE type;
    int do_check;

    G_get_window(&window);

    nrows = Rast_window_rows();
    ncols = Rast_window_cols();

    /* open the data files, input raster should be set-up already */
    if ((infd = theState->fd_old) < 0)
	G_fatal_error(_("Unable to open raster map <%s>"),
		      theState->inraster);
    if (theState->docover == TRUE) {
	if ((cinfd = theState->fd_cold) < 0)
	    G_fatal_error(_("Unable to open raster map <%s>"),
			  theState->inrcover);
    }

    if (theState->outraster != NULL) {
	if (theState->docover == TRUE)
	    type = theState->cover.type;
	else
	    type = theState->buf.type;
	outfd = Rast_open_new(theState->outraster, type);
	theState->fd_new = outfd;

    }

    if (theState->outvector) {
	if (Vect_open_new(&Out, theState->outvector, theState->z_geometry) < 0)
	    G_fatal_error(_("Unable to create vector map <%s>"),
			    theState->outvector);
	Vect_hist_command(&Out);

	fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);

	driver =
	    db_start_driver_open_database(fi->driver,
					  Vect_subst_var(fi->database, &Out));
	if (!driver)
	    G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
			  Vect_subst_var(fi->database, &Out), fi->driver);
        db_set_error_handler_driver(driver);
        
	Vect_map_add_dblink(&Out, 1, NULL, fi->table, GV_KEY_COLUMN, fi->database,
			    fi->driver);

	if (theState->docover == TRUE)
	    table = db_alloc_table(3);
	else
	    table = db_alloc_table(2);
	db_set_table_name(table, fi->table);

	column = db_get_table_column(table, 0);
	db_set_column_name(column, GV_KEY_COLUMN);
	db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);

	column = db_get_table_column(table, 1);
	db_set_column_name(column, "value");
	db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);

	if (theState->docover == TRUE) {
	    column = db_get_table_column(table, 2);
	    db_set_column_name(column, "covervalue");
	    db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
	}
	if (db_create_table(driver, table) != DB_OK)
	    G_warning(_("Cannot create new table"));

	db_begin_transaction(driver);

	Points = Vect_new_line_struct();
	Cats = Vect_new_cats_struct();
	db_init_string(&sql);
    }

    if (theState->outvector && theState->outraster)
	G_message(_("Writing raster map <%s> and vector map <%s> ..."),
		  theState->outraster, theState->outvector);
    else if (theState->outraster)
	G_message(_("Writing raster map <%s> ..."), theState->outraster);
    else if (theState->outvector)
	G_message(_("Writing vector map <%s> ..."), theState->outvector);

//.........这里部分代码省略.........

//.........这里部分代码省略.........

    /* set the pointers for multi-typed functions */
    set_func_pointers(in_type);

    /* get the window information  */
    G_get_window(&window);
    nrows = Rast_window_rows();
    ncols = Rast_window_cols();

    /* buffers for internal use */
    bndC.ns = ncols;
    bndC.sz = sizeof(CELL) * ncols;
    bndC.b[0] = G_calloc(ncols, sizeof(CELL));
    bndC.b[1] = G_calloc(ncols, sizeof(CELL));
    bndC.b[2] = G_calloc(ncols, sizeof(CELL));

    /* buffers for external use */
    bnd.ns = ncols;
    bnd.sz = ncols * bpe();
    bnd.b[0] = G_calloc(ncols, bpe());
    bnd.b[1] = G_calloc(ncols, bpe());
    bnd.b[2] = G_calloc(ncols, bpe());

    in_buf = get_buf();

    tempfile1 = G_tempfile();
    tempfile2 = G_tempfile();
    tempfile3 = G_tempfile();

    fe = open(tempfile1, O_RDWR | O_CREAT, 0666);	/* elev */
    fd = open(tempfile2, O_RDWR | O_CREAT, 0666);	/* dirn */
    fm = open(tempfile3, O_RDWR | O_CREAT, 0666);	/* problems */

    G_message(_("Reading elevation map..."));
    for (i = 0; i < nrows; i++) {
	G_percent(i, nrows, 2);
	get_row(map_id, in_buf, i);
	write(fe, in_buf, bnd.sz);
    }
    G_percent(1, 1, 1);
    Rast_close(map_id);

    /* fill single-cell holes and take a first stab at flow directions */
    G_message(_("Filling sinks..."));
    filldir(fe, fd, nrows, &bnd);

    /* determine flow directions for ambiguous cases */
    G_message(_("Determining flow directions for ambiguous cases..."));
    resolve(fd, nrows, &bndC);

    /* mark and count the sinks in each internally drained basin */
    nbasins = dopolys(fd, fm, nrows, ncols);
    if (flag1->answer) {
	/* determine the watershed for each sink */
	wtrshed(fm, fd, nrows, ncols, 4);

	/* fill all of the watersheds up to the elevation necessary for drainage */
	ppupdate(fe, fm, nrows, nbasins, &bnd, &bndC);

	/* repeat the first three steps to get the final directions */
	G_message(_("Repeat to get the final directions..."));
	filldir(fe, fd, nrows, &bnd);
	resolve(fd, nrows, &bndC);
	nbasins = dopolys(fd, fm, nrows, ncols);
    }

/* *************************************************************** */
int test_array_3d(void)
{
    int sum = 0, res = 0;
    char buff[1024];
    RASTER3D_Region region;

    N_array_3d *data1;
    N_array_3d *data11;
    N_array_3d *data2;
    N_array_3d *data22;

    N_array_3d *tmp;
    double min, max, ssum;
    int nonzero;

    /*Alloacte memory for all arrays */
    data1 =
        N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
                         FCELL_TYPE);
    N_print_array_3d_info(data1);
    data11 =
        N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
                         FCELL_TYPE);
    data2 =
        N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
                         DCELL_TYPE);
    N_print_array_3d_info(data2);
    data22 =
        N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
                         DCELL_TYPE);


    /*Fill the first arrays with data */

    res = fill_array_3d(data1);
    if (res != 0)
        G_warning("test_array_3d: error while filling array with values");
    sum += res;
    res = fill_array_3d(data2);
    if (res != 0)
        G_warning("test_array_3d: error while filling array with values");
    sum += res;

    /*Copy the data */
    N_copy_array_3d(data1, data11);
    N_copy_array_3d(data2, data22);

    /*Compare the data */
    res = compare_array_3d(data1, data11);
    if (res != 0)
        G_warning("test_array_3d: error in  N_copy_array_2d");
    sum += res;
    res = compare_array_3d(data1, data11);
    if (res != 0)
        G_warning("test_array_3d: error in  N_copy_array_2d");
    sum += res;



    /*compute statistics */
    N_calc_array_3d_stats(data1, &min, &max, &ssum, &nonzero, 0);
    G_message("FELL Min %g Max %g Sum %g  nonzero %i\n", min, max, ssum,
              nonzero);
    if (min != 0 || max != 729 || ssum != 91125 || nonzero != 1000) {
        G_warning("test_array_3d: error in  N_calc_array_3d_stats");
        sum++;
    }
    N_calc_array_3d_stats(data1, &min, &max, &ssum, &nonzero, 1);
    G_message("FELL Min %g Max %g Sum %g  nonzero %i\n", min, max, ssum,
              nonzero);
    if (min != 0 || max != 729 || ssum != 91125 || nonzero != 2744) {
        G_warning("test_array_3d: error in  N_calc_array_3d_stats");
        sum++;
    }

    N_calc_array_3d_stats(data2, &min, &max, &ssum, &nonzero, 0);
    G_message("DCELL Min %g Max %g Sum %g  nonzero %i\n", min, max, ssum,
              nonzero);
    if (min != 0 || max != 729 || ssum != 91125 || nonzero != 1000) {
        G_warning("test_array_3d: error in  N_calc_array_3d_stats");
        sum++;
    }

    N_calc_array_3d_stats(data2, &min, &max, &ssum, &nonzero, 1);
    G_message("DCELL Min %g Max %g Sum %g  nonzero %i\n", min, max, ssum,
              nonzero);
    if (min != 0 || max != 729 || ssum != 91125 || nonzero != 2744) {
        G_warning("test_array_3d: error in  N_calc_array_3d_stats");
        sum++;
    }


    /*test the array math functions */
    tmp = N_math_array_3d(data1, data2, NULL, N_ARRAY_SUM);
    N_math_array_3d(data2, data2, tmp, N_ARRAY_SUM);
    res = N_convert_array_3d_null_to_zero(tmp);
    if (res != 0)
        G_warning("test_array_3d: error in  N_convert_array_3d_null_to_zero");
    sum = res;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	Vect_close(&In);
	exit(EXIT_SUCCESS);
    }
    with_z = Vect_is_3d(&In);

    if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
	Vect_close(&In);
	G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
    }

    if (error_out->answer) {
        if (0 > Vect_open_new(&Error, error_out->answer, with_z)) {
	    Vect_close(&In);
	    G_fatal_error(_("Unable to create error vector map <%s>"), error_out->answer);
        }
    }


    Vect_copy_head_data(&In, &Out);
    Vect_hist_copy(&In, &Out);
    Vect_hist_command(&Out);

    total_input = total_output = 0;

    layer = Vect_get_field_number(&In, field_opt->answer);
    /* parse filter options */
    if (layer > 0)
	cat_list = Vect_cats_set_constraint(&In, layer, 
			      where_opt->answer, cat_opt->answer);

    if (method == DISPLACEMENT) {
	/* modifies only lines, all other features including boundaries are preserved */
	/* options where, cats, and layer are respected */
	G_message(_("Displacement..."));
	snakes_displacement(&In, &Out, thresh, alpha, beta, 1.0, 10.0,
			    iterations, cat_list, layer);
    }

    /* TODO: rearrange code below. It's really messy */
    if (method == NETWORK) {
	/* extracts lines of selected type, all other features are discarded */
	/* options where, cats, and layer are ignored */
	G_message(_("Network generalization..."));
	total_output =
	    graph_generalization(&In, &Out, mask_type, degree_thresh, 
	                         closeness_thresh, betweeness_thresh);
    }

    /* copy tables here because method == NETWORK is complete and 
     * tables for Out may be needed for parse_filter_options() below */
    if (!notab_flag->answer) {
	if (method == NETWORK)
	    copy_tables_by_cats(&In, &Out);
	else
	    Vect_copy_tables(&In, &Out, -1);
    }
    else if (where_opt->answer && method < NETWORK) {
	G_warning(_("Attributes are needed for 'where' option, copying table"));
	Vect_copy_tables(&In, &Out, -1);
    }

    /* smoothing/simplification */
    if (method < NETWORK) {
	/* modifies only lines of selected type, all other features are preserved */
	int not_modified_boundaries = 0, n_oversimplified = 0;
	struct line_pnts *APoints;  /* original Points */

//.........这里部分代码省略.........
            fprintf(stdout, _("Estimated point density: %.4g"), dens);
            fprintf(stdout, _("Estimated mean distance between points: %.4g"), dist);
        }
        else {
            fprintf(stdout, _("No points in current region"));
        }

        Vect_close(&In);
        exit(EXIT_SUCCESS);
    }

    /*----------------------------------------------------------------*/
    /* Cross-correlation begins */
    if (cross_corr_flag->answer) {
        G_debug(1, "CrossCorrelation()");
        cross = cross_correlation(&In, stepE, stepN);

        if (cross != TRUE)
            G_fatal_error(_("Cross validation didn't finish correctly"));
        else {
            G_debug(1, "Cross validation finished correctly");

            Vect_close(&In);

            G_done_msg(_("Cross validation finished for ew_step = %f and ns_step = %f"), stepE, stepN);
            exit(EXIT_SUCCESS);
        }
    }

    /* Open input ext vector */
    ext = FALSE;
    if (in_ext_opt->answer) {
        ext = TRUE;
        G_message(_("Vector map <%s> of sparse points will be interpolated"),
                  in_ext_opt->answer);

        if ((mapset = G_find_vector2(in_ext_opt->answer, "")) == NULL)
            G_fatal_error(_("Vector map <%s> not found"), in_ext_opt->answer);

        Vect_set_open_level(1);	/* WITHOUT TOPOLOGY */
        if (1 > Vect_open_old(&In_ext, in_ext_opt->answer, mapset))
            G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
                          in_opt->answer);
    }

    /* Open output map */
    /* vector output */
    if (vector && !map) {
        if (strcmp(drv, "dbf") == 0)
            G_fatal_error(_("Sorry, the <%s> driver is not compatible with "
                            "the vector output of this module. "
                            "Try with raster output or another driver."), drv);

        Vect_check_input_output_name(in_opt->answer, out_opt->answer,
                                     G_FATAL_EXIT);
        grid = FALSE;

        if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z))
            G_fatal_error(_("Unable to create vector map <%s>"),
                          out_opt->answer);

        /* Copy vector Head File */
        if (ext == FALSE) {
            Vect_copy_head_data(&In, &Out);
            Vect_hist_copy(&In, &Out);
        }

//.........这里部分代码省略.........
		if (wdepth)
		    sigmax = amax1(sigmax, sigma[k][l]);
		cchezmax = amax1(cchezmax, cchez[k][l]);
		/* saved sqrt(sinsl)*cchez to cchez array for output */
		cchez[k][l] *= sqrt(sinsl);
	    }			/* DEFined area */
	}
    }
    if (inf != NULL && smax < infmax)
	G_warning(_("Infiltration exceeds the rainfall rate everywhere! No overland flow."));

    cc = (double)mx *my;

    si0 = sisum / cc;
    vmean = vsum / cc;
    chmean = chsum / cc;

    if (inf)
	infmean = infsum / cc;

    if (wdepth)
	deltaw = 0.8 / (sigmax * vmax);	/*time step for sediment */
    deltap = 0.25 * sqrt(stepx * stepy) / vmean;	/*time step for water */

    if (deltaw > deltap)
	timec = 4.;
    else
	timec = 1.25;

    miter = (int)(timesec / (deltap * timec));	/* number of iterations = number of cells to pass */
    iterout = (int)(iterout / (deltap * timec));	/* number of cells to pass for time series output */

    fprintf(stderr, "\n");
    G_message(_("Min elevation \t= %.2f m\nMax elevation \t= %.2f m\n"), zmin,
	      zmax);
    G_message(_("Mean Source Rate (rainf. excess or sediment) \t= %f m/s or kg/m2s \n"),
	      si0);
    G_message(_("Mean flow velocity \t= %f m/s\n"), vmean);
    G_message(_("Mean Mannings \t= %f\n"), 1.0 / chmean);

    deltap = amin1(deltap, deltaw);

    G_message(_("Number of iterations \t= %d cells\n"), miter);
    G_message(_("Time step \t= %.2f s\n"), deltap);
    if (wdepth) {
	G_message(_("Sigmax \t= %f\nMax velocity \t= %f m/s\n"), sigmax,
		  vmax);
	G_message(_("Time step used \t= %.2f s\n"), deltaw);
    }
    /*    if (wdepth) deltap = 0.1; 
     *    deltap for sediment is ar. average deltap and deltaw */
    /*    if (wdepth) deltap = (deltaw+deltap)/2.; 
     *    deltap for sediment is ar. average deltap and deltaw */


    /*! For each cell (k,l) compute the length s=(v1,v2) of the path 
     *  that the particle will travel per one time step
     *  \f$ s(k,l)=v(k,l)*dt \f$, [m]=[m/s]*[s]
     *  give warning if there is a cell that will lead to path longer than 2 cells 
     *
     *  if running erosion, compute sediment transport capacity for each cell si(k,l)
     *  \f$
     * T({\bf r})=K_t({\bf r}) \bigl[\tau({\bf r})\bigr]^p
     * =K_t({\bf r}) \bigl[\rho_w\, g h({\bf r}) \sin \beta ({\bf r}) \bigr]^p
     * \f$
     * [kg/ms]=...

/* ************************************************************************* */
int main(int argc, char *argv[]) {
    struct GModule *module;
    int returnstat = 0, i;

    /* Initialize GRASS */
    G_gisinit(argv[0]);

    module = G_define_module();
    module->description
            = _("Performs unit and integration tests for the g3d library");

    /* Get parameters from user */
    set_params();

    if (G_parser(argc, argv))
        exit(EXIT_FAILURE);
    
    /* Initiate the defaults for testing */
    Rast3d_init_defaults();

    /*Run the unit tests */
    if (param.testunit->answer || param.full->answer) {
        returnstat += unit_test_coordinate_transform();
        returnstat += unit_test_put_get_value();
    }

    /*Run the integration tests */
    if (param.testint->answer || param.full->answer) {
        ;
    }

    /*Run single tests */
    if (!param.full->answer) {
        /*unit tests */
        if (!param.testunit->answer) {
            i = 0;
            if (param.unit->answers)
                while (param.unit->answers[i]) {
                    if (strcmp(param.unit->answers[i], "coord") == 0)
                        returnstat += unit_test_coordinate_transform();
                    if (strcmp(param.unit->answers[i], "putget") == 0)
                        returnstat += unit_test_put_get_value();
                    
                    i++;
                }
        }
        /*integration tests */
        if (!param.testint->answer) {
            i = 0;
            if (param.integration->answers)
                while (param.integration->answers[i]) {
                    ;
                }

        }
    }

    
    if (returnstat != 0)
        G_warning("Errors detected while testing the g3d lib");
    else
        G_message("\n-- g3d lib tests finished successfully --");

    return (returnstat);
}

int main(int argc, char *argv[])
{
    struct GModule *module;
    int Out_proj;
    int out_stat;
    int old_zone, old_proj;
    int i;
    int stat;
    char cmnd2[500];
    char proj_out[20], proj_name[50], set_name[20];
    char path[1024], buffa[1024], buffb[1024], answer[200], answer1[200];
    char answer2[200], buff[1024];
    char tmp_buff[20], *buf;

    struct Key_Value *old_proj_keys, *out_proj_keys, *in_unit_keys;
    double aa, e2;
    double f;
    FILE *FPROJ;
    int exist = 0;
    char spheroid[100];
    int j, k, sph_check;
    struct Cell_head cellhd;
    char datum[100], dat_ellps[100], dat_params[100];
    struct proj_parm *proj_parms;

    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("general"));
    G_add_keyword(_("projection"));
    module->description =
	_("Interactively reset the location's projection settings.");

    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);


    if (strcmp(G_mapset(), "PERMANENT") != 0)
	G_fatal_error(_("You must be in the PERMANENT mapset to run g.setproj"));

	/***
         * no longer necessary, table is a static struct 
	 * init_unit_table();
        ***/
    sprintf(set_name, "PERMANENT");
    G_file_name(path, "", PROJECTION_FILE, set_name);

    /* get the output projection parameters, if existing */
    /* Check for ownership here */
    stat = G__mapset_permissions(set_name);
    if (stat == 0) {
	G_fatal_error(_("PERMANENT: permission denied"));
    }
    G_get_default_window(&cellhd);
    if (-1 == G_set_window(&cellhd))
	G_fatal_error(_("Current region cannot be set"));

    if (G_get_set_window(&cellhd) == -1)
	G_fatal_error(_("Retrieving and setting region failed"));

    Out_proj = cellhd.proj;
    old_zone = cellhd.zone;
    old_proj = cellhd.proj;

    if (access(path, 0) == 0) {
	exist = 1;
	FPROJ = fopen(path, "r");
	old_proj_keys = G_fread_key_value(FPROJ);
	fclose(FPROJ);
	buf = G_find_key_value("name", old_proj_keys);
	fprintf(stderr,
		"\nWARNING: A projection file already exists for this location\n(Filename '%s')\n",
		path);
	fprintf(stderr,
		"\nThis file contains all the parameters for the location's projection:\n  %s\n",
		buf);
	fprintf(stderr,
		"\n    Overriding this information implies that the old projection parameters\n"
		"    were incorrect.  If you change the parameters, all existing data will\n"
		"    be interpreted differently by the projection software.\n%c%c%c",
		7, 7, 7);
	fprintf(stderr,
		"    GRASS will not re-project your data automatically.\n\n");

	if (!G_yes
	    (_("Would you still like to change some of the parameters?"),
	     0)) {
	    G_message(_("The projection information will not be updated"));
	    leave(SP_NOCHANGE);
	}
    }
    out_proj_keys = G_create_key_value();

    if (exist) {
	buf = G_find_key_value("zone", old_proj_keys);
	if (buf != NULL)
	    sscanf(buf, "%d", &zone);
	if (zone != old_zone) {
	    G_warning(_("Zone in default geographic region definition: %d\n"
			" is different from zone in PROJ_INFO file: %d"),
//.........这里部分代码省略.........

int
interpolate(MELEMENT rowlist[], SHORT nrows, SHORT ncols, SHORT datarows,
	    int npoints, int out_fd, int maskfd)
{
    extern CELL *cell;

    MELEMENT *Rptr;
    EW *search, *ewptr, *current_row,	/* start row for north/south search */
     *lastrow;			/* last element in search array */
    SHORT row, col;
    NEIGHBOR *nbr_head, *Nptr;
    double sum1, sum2;

    /* initialize search array and neighbors array */
    current_row = search = (EW *) G_calloc(datarows, sizeof(EW));
    lastrow = search + datarows - 1;
    nbr_head = (NEIGHBOR *) G_calloc(npoints + 1, sizeof(NEIGHBOR));
#if 0
    nbr_head->distance = maxdist;
    nbr_head->searchptr = &(nbr_head->Mptr);	/* see replace_neighbor */
#endif

    G_message(_("Interpolating raster map <%s> (%d rows)..."), output,
	      nrows);

    for (row = 0; row < nrows; row++) {	/*  loop over rows      */
	G_percent(row+1, nrows, 2);

	/* if mask occurs, read current row of the mask */
	if (mask)
	    Rast_get_c_row(maskfd, mask, row);

	/* prepare search array for next row of interpolations */
	for (ewptr = search, Rptr = rowlist; ewptr <= lastrow;
	     Rptr++, ewptr++)
	    ewptr->start = Rptr->next;	/* start at first item in row */

	for (col = 0; col < ncols; col++) {	/*  loop over columns   */

	    /* if (row != 279 && col != 209) continue; */

	    /* don't interpolate outside of the mask */
	    if (mask && mask[col] == 0) {
		cell[col] = 0;
		continue;
	    }

	    /* make a list of npoints neighboring data pts */
	    nbr_head->next = NULL;
	    if (make_neighbors_list(search, lastrow, current_row, row, col, nbr_head, npoints)) {	/* otherwise, known data value assigned */

		/* calculate value to be set for the cell from the data values
		 * of npoints closest neighboring points        */
		sum1 = sum2 = 0.0;
		Nptr = nbr_head->next;

		do {
		    sum1 += Nptr->Mptr->value / Nptr->distance;
		    sum2 += 1.0 / Nptr->distance;
		    Nptr = Nptr->next;
		} while (Nptr);	/* to end of list */

		cell[col] = (CELL) (sum1 / sum2 + .5);
		/* fprintf (stdout,"%d,%d = %d\n", col, row, cell[col]); */

		if (error_flag)	/* output interpolation error for this cell */
		    cell[col] -= mask[col];
	    }
	}			/* end of loop over columns */

	Rast_put_row(out_fd, cell, CELL_TYPE);

	/* advance current row pointer if necessary */
	if (current_row->start->y == row && current_row != lastrow)
	    ++current_row;
    }				/* end of loop over rows */

    G_free(search);

    return 0;
}

int main(int argc, char **argv)
{
    double radius;
    double fisher, david, douglas, lloyd, lloydip, morisita;
    int i, nquads, *counts;

    struct Cell_head window;
    struct GModule *module;
    struct
    {
	struct Option *input, *field, *output, *n, *r;
    } parm;
    struct
    {
	struct Flag *g;
    } flag;
    COOR *quads;

    struct Map_info Map;

    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("vector"));
    G_add_keyword(_("statistics"));
    G_add_keyword(_("point pattern"));
    module->description = _("Indices for quadrat counts of vector point lists.");

    parm.input = G_define_standard_option(G_OPT_V_INPUT);

    parm.field = G_define_standard_option(G_OPT_V_FIELD_ALL);
    
    parm.output = G_define_standard_option(G_OPT_V_OUTPUT);
    parm.output->required = NO;
    parm.output->description =
	_("Name for output quadrat centers map (number of points is written as category)");

    parm.n = G_define_option();
    parm.n->key = "nquadrats";
    parm.n->type = TYPE_INTEGER;
    parm.n->required = YES;
    parm.n->description = _("Number of quadrats");

    parm.r = G_define_option();
    parm.r->key = "radius";
    parm.r->type = TYPE_DOUBLE;
    parm.r->required = YES;
    parm.r->description = _("Quadrat radius");

    flag.g = G_define_flag();
    flag.g->key = 'g';
    flag.g->description = _("Print results in shell script style");

    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);

    sscanf(parm.n->answer, "%d", &nquads);
    sscanf(parm.r->answer, "%lf", &radius);

    G_get_window(&window);

    /* Open input */
    Vect_set_open_level(2);
    if (Vect_open_old2(&Map, parm.input->answer, "", parm.field->answer) < 0)
	G_fatal_error(_("Unable to open vector map <%s>"), parm.input->answer);

    /* Get the quadrats */
    G_message(_("Finding quadrats..."));

    quads = find_quadrats(nquads, radius, window);

    /* Get the counts per quadrat */
    G_message(_("Counting points quadrats..."));

    counts = (int *)G_malloc(nquads * (sizeof(int)));
    count_sites(quads, nquads, counts, radius, &Map,
		Vect_get_field_number(&Map, parm.field->answer));

    Vect_close(&Map);

    /* output if requested */
    if (parm.output->answer) {
	struct Map_info Out;
	struct line_pnts *Points;
	struct line_cats *Cats;

	Points = Vect_new_line_struct();
	Cats = Vect_new_cats_struct();

	if (Vect_open_new(&Out, parm.output->answer, 0) < 0)
	    G_fatal_error(_("Unable to create vector map <%s>"),
			    parm.output->answer);

	Vect_hist_command(&Out);

	for (i = 0; i < nquads; i++) {
	    Vect_reset_line(Points);
	    Vect_reset_cats(Cats);

	    Vect_append_point(Points, quads[i].x, quads[i].y, 0.0);
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    /* Global variable & function declarations */
    struct GModule *module;
    struct {
	struct Option *orig, *real, *imag;
    } opt;
    const char *Cellmap_real, *Cellmap_imag;
    const char *Cellmap_orig;
    int realfd, imagfd,  outputfd, maskfd;	/* the input and output file descriptors */
    struct Cell_head realhead, imaghead;
    DCELL *cell_real, *cell_imag;
    CELL *maskbuf;

    int i, j;			/* Loop control variables */
    int rows, cols;		/* number of rows & columns */
    long totsize;		/* Total number of data points */
    double (*data)[2];		/* Data structure containing real & complex values of FFT */

    G_gisinit(argv[0]);

    /* Set description */
    module = G_define_module();
    G_add_keyword(_("imagery"));
    G_add_keyword(_("transformation"));
    G_add_keyword(_("Fast Fourier Transform"));
    module->description =
	_("Inverse Fast Fourier Transform (IFFT) for image processing.");

    /* define options */
    opt.real = G_define_standard_option(G_OPT_R_INPUT);
    opt.real->key = "real";
    opt.real->description = _("Name of input raster map (image fft, real part)");

    opt.imag = G_define_standard_option(G_OPT_R_INPUT);
    opt.imag->key = "imaginary";
    opt.imag->description = _("Name of input raster map (image fft, imaginary part");

    opt.orig = G_define_standard_option(G_OPT_R_OUTPUT);
    opt.orig->description = _("Name for output raster map");
    
    /*call parser */
    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);

    Cellmap_real = opt.real->answer;
    Cellmap_imag = opt.imag->answer;
    Cellmap_orig = opt.orig->answer;

    /* get and compare the original window data */
    Rast_get_cellhd(Cellmap_real, "", &realhead);
    Rast_get_cellhd(Cellmap_imag, "", &imaghead);

    if (realhead.proj   != imaghead.proj   ||
	realhead.zone   != imaghead.zone   ||
	realhead.north  != imaghead.north  ||
	realhead.south  != imaghead.south  ||
	realhead.east   != imaghead.east   ||
	realhead.west   != imaghead.west   ||
	realhead.ew_res != imaghead.ew_res ||
	realhead.ns_res != imaghead.ns_res)
	G_fatal_error(_("The real and imaginary original windows did not match"));

    Rast_set_window(&realhead);	/* set the window to the whole cell map */

    /* open input raster map */
    realfd = Rast_open_old(Cellmap_real, "");
    imagfd = Rast_open_old(Cellmap_imag, "");

    /* get the rows and columns in the current window */
    rows = Rast_window_rows();
    cols = Rast_window_cols();
    totsize = rows * cols;

    /* Allocate appropriate memory for the structure containing
       the real and complex components of the FFT.  DATA[0] will
       contain the real, and DATA[1] the complex component.
     */
    data = G_malloc(rows * cols * 2 * sizeof(double));

    /* allocate the space for one row of cell map data */
    cell_real = Rast_allocate_d_buf();
    cell_imag = Rast_allocate_d_buf();
    
#define C(i, j) ((i) * cols + (j))

    /* Read in cell map values */
    G_message(_("Reading raster maps..."));
    for (i = 0; i < rows; i++) {
	Rast_get_d_row(realfd, cell_real, i);
	Rast_get_d_row(imagfd, cell_imag, i);
	for (j = 0; j < cols; j++) {
	    data[C(i, j)][0] = cell_real[j];
	    data[C(i, j)][1] = cell_imag[j];
	}
	G_percent(i+1, rows, 2);
    }

    /* close input cell maps */
    Rast_close(realfd);
//.........这里部分代码省略.........

int ram_calculate_upstream(DCELL ** distance, CELL ** dirs,
			   DCELL ** elevation, DCELL ** tmp_elevation,
			   int near)
{
    int r, c;
    int next_r, next_c;
    double easting, northing;
    double cell_easting, cell_northing;
    int i, j, k, d;
    int done;
    int counter;
    int n_inits = 0;
    double cur_dist;
    POINT *d_inits;
    double tmp_dist = 0;
    double target_elev = 0;
    size_t elevation_data_size;
    struct Cell_head window;

    Rast_get_window(&window);

    if (elevation) {
	elevation_data_size = Rast_cell_size(DCELL_TYPE);
	for (r = 0; r < nrows; ++r)
	    memcpy(tmp_elevation[r], elevation[r],
		   ncols * elevation_data_size);
    }

    for (r = 0; r < nrows; ++r)
	for (c = 0; c < ncols; ++c) {

	    for (i = 1; i < 9; ++i) {
		if (NOT_IN_REGION(i))
		    continue;	/* out of border */

		j = DIAG(i);
		next_r = NR(i);
		next_c = NC(i);
		if (dirs[next_r][next_c] == j && distance[r][c] != 0) {	/* is contributing cell */
		    distance[r][c] = -1;
		    break;
		}
	    }
	    if (distance[r][c] == 1 && dirs[r][c] > 0)
		n_inits++;
	    else if (dirs[r][c] > 0)
		distance[r][c] = -1;
	}

    d_inits = (POINT *) G_malloc(n_inits * sizeof(POINT));

    k = 0;
    for (r = 0; r < nrows; ++r)
	for (c = 0; c < ncols; ++c) {

	    if (distance[r][c] == 1) {

		distance[r][c] = 0;
		if (elevation)
		    elevation[r][c] = 0;

		d = dirs[r][c];

		if (dirs[NR(d)][NC(d)] < 0)
		    continue;

		d_inits[k].r = r;
		d_inits[k].c = c;
		d_inits[k].cur_dist = 0;


		if (elevation)
		    d_inits[k].target_elev = tmp_elevation[r][c];

		k++;
	    }
	}

    counter = n_inits = k;
    /* return 0; */
    G_message(_("Calculate upstream parameters..."));
    while (n_inits > 0) {
	k = 0;
	G_percent((counter - n_inits), counter, 10);
	for (i = 0; i < n_inits; ++i) {
	    r = d_inits[i].r;
	    c = d_inits[i].c;
	    d = dirs[r][c];
	    next_r = NR(d);
	    next_c = NC(d);
	    tmp_dist = d_inits[i].cur_dist;

	    if (elevation)
		target_elev = d_inits[i].target_elev;

	    easting = window.west + (c + 0.5) * window.ew_res;
	    northing = window.north - (r + 0.5) * window.ns_res;
	    cell_easting = window.west + (next_c + 0.5) * window.ew_res;
	    cell_northing = window.north - (next_r + 0.5) * window.ns_res;

//.........这里部分代码省略.........

/* *************************************************************** */
int test_array_2d(void)
{
    int sum = 0, res = 0;

    struct Cell_head region;
    N_array_2d *data1;
    N_array_2d *data11;
    N_array_2d *data2;
    N_array_2d *data22;
    N_array_2d *data3;
    N_array_2d *data33;
    char buff[1024];
    double min, max, ssum;
    int nonzero;

    N_array_2d *tmp;

    /*Alloacte memory for all arrays */
    data1 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, CELL_TYPE);
    N_print_array_2d_info(data1);
    data11 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, CELL_TYPE);
    data2 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, FCELL_TYPE);
    N_print_array_2d_info(data2);
    data22 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, FCELL_TYPE);
    data3 = N_alloc