← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.ros/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.ros/main.c
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       r.ros
 
5
 * AUTHOR(S):    Jianping Xu, Rutgers University, 1993
 
6
 *               Markus Neteler <neteler itc.it>
 
7
 *               Roberto Flor <flor itc.it>, Brad Douglas <rez touchofmadness.com>,
 
8
 *               Glynn Clements <glynn gclements.plus.com>, Jachym Cepicky <jachym les-ejk.cz>
 
9
 * PURPOSE:      Generates rate of spread raster map layers for wildfire modeling
 
10
 *
 
11
 * TODO: (re)move following documentation
 
12
 *
 
13
 * This raster module creates three raster map layers:
 
14
 *      1. Base (perpendicular) rate of spread (ROS);
 
15
 *      2. Maximum (forward) ROS; and
 
16
 *      3. Direction of the Maximum ROS.
 
17
 * The calculation of the two ROS values for each raster
 
18
 * cell is based on the Fortran code by Pat Andrews (1983)
 
19
 * of the Northern Forest Fire Laboratory, USDA Forest
 
20
 * Service. These three raster map layers are expected to
 
21
 * be the inputs for a seperate GRASS raster module 
 
22
 * 'r.spread'.
 
23
 * 
 
24
 * 'r.ros' can be run in two standard GRASS modes:
 
25
 * interactive and command line. For an interactive run,
 
26
 * type in 
 
27
 *      r.ros 
 
28
 * and follow the prompts; for a command line mode, 
 
29
 * type in 
 
30
 *      r.ros [-v] model=name [moisture_1h=name] 
 
31
 *                      [moisture_10h=name] 
 
32
 *                      [moisture_100h=name] 
 
33
 *                      moisture_live=name [velocity=name]
 
34
 *                      [direction=name] [slope=name] 
 
35
 *                      [aspect=name] output=name
 
36
 * where:
 
37
 *   Flag:
 
38
 *   Raster Maps:
 
39
 *      model           1-13: the standard fuel models,
 
40
 *                      all other numbers: same as barriers;
 
41
 *      moisture_1h     100*moisture_content;
 
42
 *      moisture_10h    100*moisture_content;
 
43
 *      moisture_100h   100*moisture_content;
 
44
 *      moisture_live   100*moisture_content;
 
45
 *      velocity        ft/minute;
 
46
 *      direction       degree;
 
47
 *      slope           degree;
 
48
 *      aspect          degree starting from East, anti-clockwise;
 
49
 *      output 
 
50
 *        for Base ROS  cm/min (technically not ft/min);
 
51
 *        for Max ROS   cm/min (technically not ft/min);
 
52
 *        for Direction degree.
 
53
 *
 
54
 * Note that the name given as output will be used as 
 
55
 * PREFIX for several actual raster maps. For example, if 
 
56
 * my_ros is given as the output name, 'r.ros' will 
 
57
 * actually produce three raster maps named my_ros.base, 
 
58
 * my_ros.max, and my_ros.maxdir, or even my_ros.spotdist
 
59
 * respectively. 
 
60
 *   
 
61
 * COPYRIGHT:    (C) 2000-2009 by the GRASS Development Team
 
62
 *
 
63
 *               This program is free software under the GNU General Public
 
64
 *               License (>=v2). Read the file COPYING that comes with GRASS
 
65
 *               for details.
 
66
 *
 
67
 *****************************************************************************/
 
68
 
 
69
#include <stdlib.h>
 
70
#include <stdio.h>
 
71
#include <math.h>
 
72
#include <grass/gis.h>
 
73
#include <grass/raster.h>
 
74
#include <grass/glocale.h>
 
75
#include "local_proto.h"
 
76
 
 
77
 
 
78
#define DATA(map, r, c)         (map)[(r) * ncols + (c)]
 
79
 
 
80
/*measurements of the 13 fuel models, input of Rothermel equation (1972) */
 
81
float WO[4][14] =
 
82
    { {0, 0.034, 0.092, 0.138, 0.230, 0.046, 0.069, 0.052, 0.069, 0.134,
 
83
       0.138, 0.069, 0.184, 0.322},
 
84
{0, 0, 0.046, 0, 0.184, 0.023, 0.115, 0.086, 0.046, 0.019, 0.092, 0.207,
 
85
 0.644, 1.058},
 
86
{0, 0, 0.023, 0, 0.092, 0, 0.092, 0.069, 0.115, 0.007, 0.230, 0.253, 0.759,
 
87
 1.288},
 
88
{0, 0, 0.023, 0, 0.230, 0.092, 0, 0.017, 0, 0, 0.092, 0, 0}
 
89
};
 
90
 
 
91
                                            /*ovendry fuel loading, lb./ft.^2 */
 
92
float DELTA[] = { 0, 1.0, 1.0, 2.5, 6.0, 2.0, 2.5, 2.5,
 
93
    0.2, 0.2, 1.0, 1.0, 2.3, 3.0
 
94
};                              /*fuel depth, ft. */
 
95
float SIGMA[4][14] =
 
96
    { {0, 3500, 3000, 1500, 2000, 2000, 1750, 1750, 2000, 2500, 2000, 1500,
 
97
       1500, 1500},
 
98
{0, 0, 109, 0, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109},
 
99
{0, 0, 30, 0, 30, 0, 30, 30, 30, 30, 30, 30, 30, 30},
 
100
{0, 0, 1500, 0, 1500, 1500, 0, 1500, 0, 0, 1500, 0, 0, 0}
 
101
};
 
102
 
 
103
                         /*fuel particale surface-area-to-volume ratio, 1/ft. */
 
104
float MX[] = { 0, 0.12, 0.15, 0.25, 0.20, 0.20, 0.25, 0.40,
 
105
    0.30, 0.25, 0.25, 0.15, 0.20, 0.25
 
106
};                              /*moisture content of extinction */
 
107
 
 
108
CELL *map_elev;                 /*full array for elevation map layer (for spotting) */
 
109
int nrows, ncols;
 
110
struct Cell_head window;
 
111
 
 
112
 
 
113
int main(int argc, char *argv[])
 
114
{
 
115
 
 
116
        /***input of Rothermel equation (1972)***/
 
117
    float h = 8000.0,           /*heat of combustion, BTU/lb. */
 
118
        rhop = 32,              /*ovendry fuel density, lb./ft.^3 */
 
119
        ST = 0.0555;            /*fuel total mineral content, lb. minerals/lb. ovendry */
 
120
 
 
121
    float sigma[14];
 
122
 
 
123
        /***derived parameters of Rothermel equation (1972)***/
 
124
    float R,                    /*rate of spread, ft./min.
 
125
                                   R = IR*xi*(1+phiw+phis)/(rhob*epsilon*Qig)  */
 
126
      IR,                       /*reaction intensity, BTU/ft.^2/min.
 
127
                                   IR = gamma*wn*h*etaM*etas  */
 
128
      gamma,                    /*optimum reation velosity, 1/min.
 
129
                                   gamma = gammamax*(beta/betaop)^A*
 
130
                                   exp(A(1-beta/betaop))  */
 
131
      gammamax,                 /*maximum reation velosity, 1/min.
 
132
                                   gammamax = sigma^1.5/(495+0.0594*sigma^1.5)  */
 
133
      betaop,                   /*optimum packing ratio 
 
134
                                   betaop = 3.348/(sigma^0.8189)  */
 
135
      A,                        /*A = 1/(4.77*sigma^0.1-7.27)  */
 
136
      etas = 0.174 / pow(0.01, 0.19),   /*mineral damping coefficient */
 
137
        xi,                     /*propagating flux ratio,
 
138
                                   xi = exp((0.792+0.681*sigma^0.5)(beta+0.1))/
 
139
                                   (192+0.2595*sigma)  */
 
140
        phiw,                   /*wind coefficient,  phiw = C*U^B/(beta/betaop)^E */
 
141
        C,                      /*C = 7.47/exp(0.133*sigma^0.55)  */
 
142
        B,                      /*B = 0.02526*sigma^.054   */
 
143
        E,                      /*E = 0.715/exp(0.000359*sigma)  */
 
144
        phis,                   /*slope coefficient,phis = 5.275/beta^0.3*(tan(theta)^2) */
 
145
        rhob,                   /*ovendry bulk density, lb./ft.^3, rohb = wo/delta */
 
146
        epsilon[4][14],         /*effective heating number, epsilon = 1/exp(138/sigma) */
 
147
        Qig[14],                /*heat of preignition, BTU/lb.  Qig = 250+1116*Mf */
 
148
        beta;                   /*packing ratio,  beta = rhob/rhop  */
 
149
 
 
150
        /***intermediate variables***/
 
151
    float R0,                   /*base ROS (w/out wind/slope) */
 
152
      Rdir, sin_fac, cos_fac, Ffactor_all[4][14],       /*in all fuel subclasses by sigma/WO */
 
153
      Ffactor_in_dead[3][14],   /*in dead fuel subclasses by sigma/WO */
 
154
      Ffactor_dead[14],         /*dead fuel weight by sigma/WO */
 
155
      Ffactor_live[14],         /*live fuel weight by sigma/WO */
 
156
      Gfactor_in_dead[3][14],   /*in dead fuel by the 6 classes */
 
157
      G1, G2, G3, G4, G5, wo_dead[14],  /*dead fuel total load */
 
158
      wn_dead,                  /*net dead fuel total load */
 
159
      wn_live,                  /*net live fuel (total) load */
 
160
      class_sum, moisture[4],   /*moistures of 1-h,10-h,100-h,live fuels */
 
161
      Mf_dead,                  /*total moisture of dead fuels */
 
162
      etaM_dead,                /*dead fuel misture damping coefficent */
 
163
      etaM_live,                /*live fuel misture damping coefficent */
 
164
      xmext,                    /*live fuel moisture of extinction */
 
165
      phi_ws,                   /*wind and slope conbined coefficient */
 
166
      wmfd, fdmois, fined, finel;
 
167
 
 
168
    /*other local variables */
 
169
    int col, row,
 
170
        spotting,
 
171
        model, class,
 
172
        fuel_fd = 0,
 
173
        mois_1h_fd = 0, mois_10h_fd = 0, mois_100h_fd = 0, mois_live_fd = 0,
 
174
        vel_fd = 0, dir_fd = 0,
 
175
        elev_fd = 0, slope_fd = 0, aspect_fd = 0,
 
176
        base_fd = 0, max_fd = 0, maxdir_fd = 0, spotdist_fd = 0;
 
177
 
 
178
    char *name_base, *name_max, *name_maxdir, *name_spotdist;
 
179
 
 
180
    CELL *fuel,                 /*cell buffer for fuel model map layer */
 
181
     *mois_1h,                  /*cell buffer for 1-hour fuel moisture map layer */
 
182
     *mois_10h,                 /*cell buffer for 10-hour fuel moisture map layer */
 
183
     *mois_100h,                /*cell buffer for 100-hour fuel moisture map layer */
 
184
     *mois_live,                /*cell buffer for live fuel moisture map layer */
 
185
     *vel,                      /*cell buffer for wind velocity map layer */
 
186
     *dir,                      /*cell buffer for wind direction map layer */
 
187
     *elev = NULL,              /*cell buffer for elevation map layer (for spotting) */
 
188
        *slope,                 /*cell buffer for slope map layer */
 
189
        *aspect,                /*cell buffer for aspect map layer */
 
190
        *base,                  /*cell buffer for base ROS map layer */
 
191
        *max,                   /*cell buffer for max ROS map layer */
 
192
        *maxdir,                /*cell buffer for max ROS direction map layer */
 
193
        *spotdist = NULL;       /*cell buffer for max spotting distance map layer */
 
194
 
 
195
    extern struct Cell_head window;
 
196
 
 
197
    struct
 
198
    {
 
199
        struct Option *model,
 
200
            *mois_1h, *mois_10h, *mois_100h, *mois_live,
 
201
            *vel, *dir, *elev, *slope, *aspect, *base, *max, *maxdir, *spotdist;
 
202
    } parm;
 
203
 
 
204
    /* please, remove before GRASS 7 released */
 
205
    struct GModule *module;
 
206
 
 
207
    /* initialize access to database and create temporary files */
 
208
    G_gisinit(argv[0]);
 
209
 
 
210
    /* Set description */
 
211
    module = G_define_module();
 
212
    G_add_keyword(_("raster"));
 
213
    G_add_keyword(_("fire"));
 
214
    G_add_keyword(_("spread"));
 
215
    G_add_keyword(_("rate of spread"));
 
216
    G_add_keyword(_("hazard"));
 
217
    module->label = _("Generates rate of spread raster maps.");
 
218
    module->description =
 
219
        _("Generates three, or four raster map layers showing the base "
 
220
          "(perpendicular) rate of spread (ROS), the maximum (forward) ROS, "
 
221
          "the direction of the maximum ROS, and optionally the "
 
222
          "maximum potential spotting distance for fire spread simulation.");
 
223
 
 
224
    parm.model = G_define_standard_option(G_OPT_R_INPUT);
 
225
    parm.model->key = "model";
 
226
    parm.model->label = _("Raster map containing fuel models");
 
227
    parm.model->description =
 
228
        _("Name of an "
 
229
          "existing raster map layer in the user's current mapset search path containing "
 
230
          "the standard fuel models defined by the USDA Forest Service. Valid values "
 
231
          "are 1-13; other numbers are recognized as barriers by r.ros.");
 
232
 
 
233
    parm.mois_1h = G_define_standard_option(G_OPT_R_INPUT);
 
234
    parm.mois_1h->key = "moisture_1h";
 
235
    parm.mois_1h->required = NO;
 
236
    parm.mois_1h->label =
 
237
        _("Raster map containing the 1-hour fuel moisture (%)");
 
238
        parm.mois_1h->description =
 
239
        _("Name of an existing raster map layer in "
 
240
          "the user's current mapset search path containing the 1-hour (<.25\") "
 
241
          "fuel moisture (percentage content multiplied by 100).");
 
242
 
 
243
    parm.mois_10h = G_define_standard_option(G_OPT_R_INPUT);
 
244
    parm.mois_10h->key = "moisture_10h";
 
245
    parm.mois_10h->required = NO;
 
246
    parm.mois_10h->label =
 
247
        _("Raster map containing the 10-hour fuel moisture (%)");
 
248
    parm.mois_10h->description =
 
249
        _("Name of an existing raster map layer in the "
 
250
          "user's current mapset search path containing the 10-hour (.25-1\") fuel "
 
251
          "moisture (percentage content multiplied by 100).");
 
252
 
 
253
    parm.mois_100h = G_define_standard_option(G_OPT_R_INPUT);
 
254
    parm.mois_100h->key = "moisture_100h";
 
255
    parm.mois_100h->required = NO;
 
256
    parm.mois_100h->label =
 
257
        _("Raster map containing the 100-hour fuel moisture (%)");
 
258
    parm.mois_100h->description =
 
259
        _("Name of an existing raster map layer in the "
 
260
          "user's current mapset search path containing the 100-hour (1-3\") fuel moisture "
 
261
          "(percentage content multiplied by 100).");
 
262
 
 
263
    parm.mois_live = G_define_standard_option(G_OPT_R_INPUT);
 
264
    parm.mois_live->key = "moisture_live";
 
265
    parm.mois_live->label =
 
266
        _("Raster map containing live fuel moisture (%)");
 
267
    parm.mois_live->description =
 
268
        _("Name of an existing raster map layer in the "
 
269
          "user's current mapset search path containing live (herbaceous) fuel "
 
270
          "moisture (percentage content multiplied by 100).");
 
271
 
 
272
    parm.vel = G_define_standard_option(G_OPT_R_INPUT);
 
273
    parm.vel->key = "velocity";
 
274
    parm.vel->required = NO;
 
275
    parm.vel->label =
 
276
        _("Raster map containing midflame wind velocities (ft/min)");
 
277
    parm.vel->description =
 
278
        _("Name of an existing raster map layer in the user's "
 
279
          "current mapset search path containing wind velocities at half of the average "
 
280
          "flame height (feet/minute).");
 
281
 
 
282
    parm.dir = G_define_standard_option(G_OPT_R_INPUT);
 
283
    parm.dir->key = "direction";
 
284
    parm.dir->required = NO;
 
285
    parm.dir->label =
 
286
        _("Name of raster map containing wind directions (degree)");
 
287
    parm.dir->description =
 
288
        _("Name of an existing raster map "
 
289
          "layer in the user's current mapset search path containing wind direction, "
 
290
          "clockwise from north (degree).");
 
291
 
 
292
    parm.slope = G_define_standard_option(G_OPT_R_INPUT);
 
293
    parm.slope->key = "slope";
 
294
    parm.slope->required = NO;
 
295
    parm.slope->label =
 
296
        _("Name of raster map containing slope (degree)");
 
297
    parm.slope->description =
 
298
        _("Name of an existing raster map layer "
 
299
          "in the user's current mapset search path containing "
 
300
          "topographic slope (degree).");
 
301
 
 
302
    parm.aspect = G_define_standard_option(G_OPT_R_INPUT);
 
303
    parm.aspect->key = "aspect";
 
304
    parm.aspect->required = NO;
 
305
    parm.aspect->label =
 
306
        _("Raster map containing aspect (degree, CCW from E)");
 
307
    parm.aspect->description =
 
308
        _("Name of an existing "
 
309
          "raster map layer in the user's current mapset search path containing "
 
310
          "topographic aspect, counterclockwise from east (GRASS convention) "
 
311
          "in degrees.");
 
312
 
 
313
    parm.elev = G_define_standard_option(G_OPT_R_ELEV);
 
314
    parm.elev->required = NO;
 
315
    parm.elev->label =
 
316
        _("Raster map containing elevation (m, required for spotting)");
 
317
    parm.elev->description =
 
318
        _("Name of an existing raster map "
 
319
          "layer in the user's current mapset search path containing elevation (meters). "
 
320
          "Option is required from spotting distance computation "
 
321
          "(when spotting_distance option is provided)");
 
322
 
 
323
    parm.base = G_define_standard_option(G_OPT_R_OUTPUT);
 
324
    parm.base->key = "base_ros";
 
325
    parm.base->required = YES;
 
326
    parm.base->label =
 
327
        _("Output raster map containing base ROS (cm/min)");
 
328
    parm.base->description =
 
329
        _("Base (perpendicular) rate of spread (ROS)");
 
330
 
 
331
    parm.max = G_define_standard_option(G_OPT_R_OUTPUT);
 
332
    parm.max->key = "max_ros";
 
333
    parm.max->required = YES;
 
334
    parm.max->label =
 
335
        _("Output raster map containing maximal ROS (cm/min)");
 
336
    parm.max->description =
 
337
        _("The maximum (forward) rate of spread (ROS)");
 
338
 
 
339
    parm.maxdir = G_define_standard_option(G_OPT_R_OUTPUT);
 
340
    parm.maxdir->key = "direction_ros";
 
341
    parm.maxdir->required = YES;
 
342
    parm.maxdir->label =
 
343
        _("Output raster map containing directions of maximal ROS (degree)");
 
344
    parm.maxdir->description =
 
345
        _("The direction of the maximal (forward) rate of spread (ROS)");
 
346
 
 
347
    parm.spotdist = G_define_standard_option(G_OPT_R_OUTPUT);
 
348
    parm.spotdist->key = "spotting_distance";
 
349
    parm.spotdist->required = NO;
 
350
    parm.spotdist->label =
 
351
        _("Output raster map containing maximal spotting distance (m)");
 
352
    parm.spotdist->description =
 
353
        _("The maximal potential spotting distance"
 
354
          " (requires elevation raster map to be provided).");
 
355
 
 
356
    /*   Parse command line */
 
357
    if (G_parser(argc, argv))
 
358
        exit(EXIT_FAILURE);
 
359
 
 
360
    if (parm.spotdist->answer)
 
361
        spotting = 1;
 
362
    else
 
363
        spotting = 0;
 
364
 
 
365
    /*  Check if input layers exists in data base  */
 
366
    if (G_find_raster2(parm.model->answer, "") == NULL)
 
367
        G_fatal_error(_("Raster map <%s> not found"), parm.model->answer);
 
368
 
 
369
    if (!
 
370
        (parm.mois_1h->answer || parm.mois_10h->answer ||
 
371
         parm.mois_100h->answer)) {
 
372
        G_fatal_error(_("No dead fuel moisture is given. "
 
373
                        "At least one of the 1-h, 10-h, 100-h moisture layers is required."));
 
374
    }
 
375
 
 
376
    if (parm.mois_1h->answer) {
 
377
        if (G_find_raster2(parm.mois_1h->answer, "") == NULL)
 
378
            G_fatal_error(_("Raster map <%s> not found"),
 
379
                          parm.mois_1h->answer);
 
380
    }
 
381
    if (parm.mois_10h->answer) {
 
382
        if (G_find_raster2(parm.mois_10h->answer, "") == NULL)
 
383
            G_fatal_error(_("Raster map <%s> not found"),
 
384
                          parm.mois_10h->answer);
 
385
    }
 
386
    if (parm.mois_100h->answer) {
 
387
        if (G_find_raster2(parm.mois_100h->answer, "") == NULL)
 
388
            G_fatal_error(_("Raster map <%s> not found"),
 
389
                          parm.mois_100h->answer);
 
390
    }
 
391
 
 
392
    if (G_find_raster2(parm.mois_live->answer, "") == NULL)
 
393
        G_fatal_error(_("Raster map <%s> not found"), parm.mois_live->answer);
 
394
 
 
395
    if (parm.vel->answer && !(parm.dir->answer)) {
 
396
        G_fatal_error(_("A wind direction layer should be "
 
397
                        "given if the wind velocity layer <%s> has been given"),
 
398
                      parm.vel->answer);
 
399
    }
 
400
    if (!(parm.vel->answer) && parm.dir->answer) {
 
401
        G_fatal_error(_("A wind velocity layer should be given "
 
402
                        "if the wind direction layer <%s> has been given"),
 
403
                      parm.dir->answer);
 
404
    }
 
405
    if (parm.vel->answer) {
 
406
        if (G_find_raster2(parm.vel->answer, "") == NULL)
 
407
            G_fatal_error(_("Raster map <%s> not found"), parm.vel->answer);
 
408
    }
 
409
    if (parm.dir->answer) {
 
410
        if (G_find_raster2(parm.dir->answer, "") == NULL)
 
411
            G_fatal_error(_("Raster map <%s> not found"), parm.dir->answer);
 
412
    }
 
413
 
 
414
    if (parm.slope->answer && !(parm.aspect->answer)) {
 
415
        G_fatal_error(_("An aspect layer should be given "
 
416
                        "if the slope layer <%s> has been given"),
 
417
                      parm.slope->answer);
 
418
    }
 
419
    if (!(parm.slope->answer) && parm.aspect->answer) {
 
420
        G_fatal_error(_("A slope layer should be given "
 
421
                        "if the aspect layer <%s> has been given"),
 
422
                      parm.aspect->answer);
 
423
    }
 
424
    if (parm.slope->answer) {
 
425
        if (G_find_raster2(parm.slope->answer, "") == NULL)
 
426
            G_fatal_error(_("Raster map <%s> not found"), parm.slope->answer);
 
427
    }
 
428
    if (parm.aspect->answer) {
 
429
        if (G_find_raster2(parm.aspect->answer, "") == NULL)
 
430
            G_fatal_error(_("Raster map <%s> not found"),
 
431
                          parm.aspect->answer);
 
432
    }
 
433
 
 
434
    if (spotting) {
 
435
        if (!(parm.elev->answer)) {
 
436
            G_fatal_error(_("An elevation layer should be given "
 
437
                            "if considering spotting"));
 
438
        }
 
439
        else {
 
440
            if (G_find_raster2(parm.elev->answer, "") == NULL)
 
441
                G_fatal_error(_("Raster map <%s> not found"),
 
442
                              parm.elev->answer);
 
443
        }
 
444
    }
 
445
 
 
446
    /*assign names of the three output ROS layers */
 
447
    name_base = parm.base->answer;
 
448
    name_max = parm.max->answer;
 
449
    name_maxdir = parm.maxdir->answer;
 
450
 
 
451
    /*assign a name to output SPOTTING distance layer */
 
452
    if (spotting) {
 
453
        name_spotdist = parm.spotdist->answer;
 
454
    }
 
455
 
 
456
    /*  Get database window parameters  */
 
457
    G_get_window(&window);
 
458
 
 
459
    /*  find number of rows and columns in window    */
 
460
    nrows = Rast_window_rows();
 
461
    ncols = Rast_window_cols();
 
462
 
 
463
    fuel = Rast_allocate_c_buf();
 
464
    mois_1h = Rast_allocate_c_buf();
 
465
    mois_10h = Rast_allocate_c_buf();
 
466
    mois_100h = Rast_allocate_c_buf();
 
467
    mois_live = Rast_allocate_c_buf();
 
468
    vel = Rast_allocate_c_buf();
 
469
    dir = Rast_allocate_c_buf();
 
470
    slope = Rast_allocate_c_buf();
 
471
    aspect = Rast_allocate_c_buf();
 
472
    base = Rast_allocate_c_buf();
 
473
    max = Rast_allocate_c_buf();
 
474
    maxdir = Rast_allocate_c_buf();
 
475
    if (spotting) {
 
476
        spotdist = Rast_allocate_c_buf();
 
477
        elev = Rast_allocate_c_buf();
 
478
        map_elev = (CELL *) G_calloc(nrows * ncols, sizeof(CELL));
 
479
    }
 
480
 
 
481
    /*  Open input cell layers for reading  */
 
482
 
 
483
    fuel_fd = Rast_open_old(parm.model->answer, "");
 
484
 
 
485
    if (parm.mois_1h->answer)
 
486
        mois_1h_fd = Rast_open_old(parm.mois_1h->answer, "");
 
487
 
 
488
    if (parm.mois_10h->answer)
 
489
        mois_10h_fd = Rast_open_old(parm.mois_10h->answer,"");
 
490
 
 
491
    if (parm.mois_100h->answer)
 
492
        mois_100h_fd = Rast_open_old(parm.mois_100h->answer, "");
 
493
 
 
494
    mois_live_fd = Rast_open_old(parm.mois_live->answer, "");
 
495
 
 
496
    if (parm.vel->answer)
 
497
        vel_fd = Rast_open_old(parm.vel->answer, "");
 
498
 
 
499
    if (parm.dir->answer)
 
500
        dir_fd = Rast_open_old(parm.dir->answer, "");
 
501
 
 
502
    if (parm.slope->answer)
 
503
        slope_fd = Rast_open_old(parm.slope->answer, "");
 
504
 
 
505
    if (parm.aspect->answer)
 
506
        aspect_fd = Rast_open_old(parm.aspect->answer, "");
 
507
 
 
508
    if (spotting)
 
509
        elev_fd = Rast_open_old(parm.elev->answer, "");
 
510
 
 
511
    base_fd = Rast_open_c_new(name_base);
 
512
    max_fd = Rast_open_c_new(name_max);
 
513
    maxdir_fd = Rast_open_c_new(name_maxdir);
 
514
    if (spotting)
 
515
        spotdist_fd = Rast_open_c_new(name_spotdist);
 
516
 
 
517
    /*compute weights, combined wo, and combined sigma */
 
518
    /*wo[model] -- simple sum of WO[class][model] by all fuel subCLASS */
 
519
    /*sigma[model] -- weighted sum of SIGMA[class][model] by all fuel subCLASS       *epsilon[class][model] */
 
520
    for (model = 1; model <= 13; model++) {
 
521
        class_sum = 0.0;
 
522
        wo_dead[model] = 0.0;
 
523
        sigma[model] = 0.0;
 
524
        for (class = 0; class <= 3; class++) {
 
525
            class_sum = class_sum + WO[class][model] * SIGMA[class][model];
 
526
            if (SIGMA[class][model] > 0.0) {
 
527
                epsilon[class][model] = exp(-138.0 / SIGMA[class][model]);
 
528
            }
 
529
            else {
 
530
                epsilon[class][model] = 0.0;
 
531
            }
 
532
        }
 
533
        for (class = 0; class <= 3; class++) {
 
534
            Ffactor_all[class][model] =
 
535
                WO[class][model] * SIGMA[class][model] / class_sum;
 
536
            sigma[model] =
 
537
                sigma[model] +
 
538
                SIGMA[class][model] * Ffactor_all[class][model];
 
539
        }
 
540
        class_sum = 0.0;
 
541
        for (class = 0; class <= 2; class++) {
 
542
            wo_dead[model] = wo_dead[model] + WO[class][model];
 
543
            class_sum = class_sum + WO[class][model] * SIGMA[class][model];
 
544
        }
 
545
        for (class = 0; class <= 2; class++) {
 
546
            Ffactor_in_dead[class][model] =
 
547
                WO[class][model] * SIGMA[class][model] / class_sum;
 
548
        }
 
549
 
 
550
        /* compute G factor for each of the 6 subclasses */
 
551
        G1 = 0.0;
 
552
        G2 = 0.0;
 
553
        G3 = 0.0;
 
554
        G4 = 0.0;
 
555
        G5 = 0.0;
 
556
        for (class = 0; class <= 2; class++) {
 
557
            if (SIGMA[class][model] >= 1200)
 
558
                G1 = G1 + Ffactor_in_dead[class][model];
 
559
            if (SIGMA[class][model] < 1200 && SIGMA[class][model] >= 192)
 
560
                G2 = G2 + Ffactor_in_dead[class][model];
 
561
            if (SIGMA[class][model] < 192 && SIGMA[class][model] >= 96)
 
562
                G3 = G3 + Ffactor_in_dead[class][model];
 
563
            if (SIGMA[class][model] < 96 && SIGMA[class][model] >= 48)
 
564
                G4 = G4 + Ffactor_in_dead[class][model];
 
565
            if (SIGMA[class][model] < 48 && SIGMA[class][model] >= 16)
 
566
                G5 = G5 + Ffactor_in_dead[class][model];
 
567
        }
 
568
        for (class = 0; class <= 2; class++) {
 
569
            if (SIGMA[class][model] >= 1200)
 
570
                Gfactor_in_dead[class][model] = G1;
 
571
            if (SIGMA[class][model] < 1200 && SIGMA[class][model] >= 192)
 
572
                Gfactor_in_dead[class][model] = G2;
 
573
            if (SIGMA[class][model] < 192 && SIGMA[class][model] >= 96)
 
574
                Gfactor_in_dead[class][model] = G3;
 
575
            if (SIGMA[class][model] < 96 && SIGMA[class][model] >= 48)
 
576
                Gfactor_in_dead[class][model] = G4;
 
577
            if (SIGMA[class][model] < 48 && SIGMA[class][model] >= 16)
 
578
                Gfactor_in_dead[class][model] = G5;
 
579
            if (SIGMA[class][model] < 16)
 
580
                Gfactor_in_dead[class][model] = 0.0;
 
581
        }
 
582
 
 
583
        Ffactor_dead[model] =
 
584
            class_sum / (class_sum + WO[3][model] * SIGMA[3][model]);
 
585
        Ffactor_live[model] = 1 - Ffactor_dead[model];
 
586
    }
 
587
 
 
588
    /*if considering spotting, read elevation map into an array */
 
589
    if (spotting)
 
590
        for (row = 0; row < nrows; row++) {
 
591
            Rast_get_c_row(elev_fd, elev, row);
 
592
            for (col = 0; col < ncols; col++)
 
593
                DATA(map_elev, row, col) = elev[col];
 
594
        }
 
595
 
 
596
    /*major computation: compute ROSs one cell a time */
 
597
    for (row = 0; row < nrows; row++) {
 
598
        G_percent(row, nrows, 2);
 
599
        Rast_get_c_row(fuel_fd, fuel, row);
 
600
        if (parm.mois_1h->answer)
 
601
            Rast_get_c_row(mois_1h_fd, mois_1h, row);
 
602
        if (parm.mois_10h->answer)
 
603
            Rast_get_c_row(mois_10h_fd, mois_10h, row);
 
604
        if (parm.mois_100h->answer)
 
605
            Rast_get_c_row(mois_100h_fd, mois_100h, row);
 
606
        Rast_get_c_row(mois_live_fd, mois_live, row);
 
607
        if (parm.vel->answer)
 
608
            Rast_get_c_row(vel_fd, vel, row);
 
609
        if (parm.dir->answer)
 
610
            Rast_get_c_row(dir_fd, dir, row);
 
611
        if (parm.slope->answer)
 
612
            Rast_get_c_row(slope_fd, slope, row);
 
613
        if (parm.aspect->answer)
 
614
            Rast_get_c_row(aspect_fd, aspect, row);
 
615
 
 
616
        /*initialize cell buffers for output map layers */
 
617
        for (col = 0; col < ncols; col++) {
 
618
            base[col] = max[col] = maxdir[col] = 0;
 
619
            if (spotting)
 
620
                spotdist[col] = 0;
 
621
        }
 
622
 
 
623
        for (col = 0; col < ncols; col++) {
 
624
            /*check if a fuel is within the 13 models,
 
625
             *if not, no processing; useful when no data presents*/
 
626
            if (fuel[col] < 1 || fuel[col] > 13)
 
627
                continue;
 
628
            if (parm.mois_1h->answer)
 
629
                moisture[0] = 0.01 * mois_1h[col];
 
630
            if (parm.mois_10h->answer)
 
631
                moisture[1] = 0.01 * mois_10h[col];
 
632
            if (parm.mois_100h->answer)
 
633
                moisture[2] = 0.01 * mois_100h[col];
 
634
            moisture[3] = 0.01 * mois_live[col];
 
635
            if (parm.aspect->answer)
 
636
                aspect[col] = (630 - aspect[col]) % 360;
 
637
 
 
638
            /* assign some dead fuel moisture if not completely 
 
639
             *     given based on emperical relationship*/
 
640
            if (!(parm.mois_10h->answer || parm.mois_100h->answer)) {
 
641
                moisture[1] = moisture[0] + 0.01;
 
642
                moisture[2] = moisture[0] + 0.02;
 
643
            }
 
644
            if (!(parm.mois_1h->answer || parm.mois_100h->answer)) {
 
645
                moisture[0] = moisture[1] - 0.01;
 
646
                moisture[2] = moisture[1] + 0.01;
 
647
            }
 
648
            if (!(parm.mois_1h->answer || parm.mois_10h->answer)) {
 
649
                moisture[0] = moisture[2] - 0.02;
 
650
                moisture[1] = moisture[2] - 0.01;
 
651
            }
 
652
            if (!(parm.mois_1h->answer) && parm.mois_10h->answer &&
 
653
                parm.mois_100h->answer)
 
654
                moisture[0] = moisture[1] - 0.01;
 
655
            if (!(parm.mois_10h->answer) && parm.mois_1h->answer &&
 
656
                parm.mois_100h->answer)
 
657
                moisture[1] = moisture[0] + 0.01;
 
658
            if (!(parm.mois_100h->answer) && parm.mois_1h->answer &&
 
659
                parm.mois_10h->answer)
 
660
                moisture[2] = moisture[1] + 0.01;
 
661
 
 
662
            /*compute xmext, moisture of extinction of live fuels */
 
663
            wmfd = 0.0;
 
664
            fined = 0.0;
 
665
            if (SIGMA[3][fuel[col]] > 0.0) {
 
666
                for (class = 0; class <= 2; class++) {
 
667
                    if (SIGMA[class][fuel[col]] == 0.0)
 
668
                        continue;
 
669
                    fined =
 
670
                        fined +
 
671
                        WO[class][fuel[col]] * exp(-138.0 /
 
672
                                                   SIGMA[class][fuel[col]]);
 
673
                    wmfd =
 
674
                        wmfd +
 
675
                        WO[class][fuel[col]] * exp(-138.0 /
 
676
                                                   SIGMA[class][fuel[col]]) *
 
677
                        moisture[class];
 
678
                }
 
679
                fdmois = wmfd / fined;
 
680
                finel = WO[3][fuel[col]] * exp(-500.0 / SIGMA[3][fuel[col]]);
 
681
                xmext =
 
682
                    2.9 * (fined / finel) * (1 - fdmois / MX[fuel[col]]) -
 
683
                    0.226;
 
684
            }
 
685
            else
 
686
                xmext = MX[fuel[col]];
 
687
            if (xmext < MX[fuel[col]])
 
688
                xmext = MX[fuel[col]];
 
689
 
 
690
            /*compute other intermediate values */
 
691
            Mf_dead = 0.0;
 
692
            wn_dead = 0.0;
 
693
            class_sum = 0.0;
 
694
            for (class = 0; class <= 2; class++) {
 
695
                Mf_dead =
 
696
                    Mf_dead +
 
697
                    moisture[class] * Ffactor_in_dead[class][fuel[col]];
 
698
                wn_dead =
 
699
                    wn_dead +
 
700
                    WO[class][fuel[col]] * Gfactor_in_dead[class][fuel[col]] *
 
701
                    (1 - ST);
 
702
                Qig[class] = 250 + 1116 * moisture[class];
 
703
                class_sum =
 
704
                    class_sum +
 
705
                    Ffactor_all[class][fuel[col]] *
 
706
                    epsilon[class][fuel[col]] * Qig[class];
 
707
            }
 
708
            etaM_dead =
 
709
                1.0 - 2.59 * (Mf_dead / MX[fuel[col]]) +
 
710
                5.11 * (Mf_dead / MX[fuel[col]]) * (Mf_dead / MX[fuel[col]]) -
 
711
                3.52 * (Mf_dead / MX[fuel[col]]) * (Mf_dead / MX[fuel[col]]) *
 
712
                (Mf_dead / MX[fuel[col]]);
 
713
            if (Mf_dead >= MX[fuel[col]])
 
714
                etaM_dead = 0.0;
 
715
            etaM_live =
 
716
                1.0 - 2.59 * (moisture[3] / xmext) +
 
717
                5.11 * (moisture[3] / xmext) * (moisture[3] / xmext) -
 
718
                3.52 * (moisture[3] / xmext) * (moisture[3] / xmext) *
 
719
                (moisture[3] / xmext);
 
720
            if (moisture[3] >= xmext)
 
721
                etaM_live = 0.0;
 
722
            wn_live = WO[3][fuel[col]] * (1 - ST);
 
723
            Qig[3] = 250 + 1116 * moisture[3];
 
724
            class_sum =
 
725
                class_sum +
 
726
                Ffactor_all[3][fuel[col]] * epsilon[3][fuel[col]] * Qig[3];
 
727
 
 
728
            /*final computations */
 
729
            rhob = (wo_dead[fuel[col]] + WO[3][fuel[col]]) / DELTA[fuel[col]];
 
730
            beta = rhob / rhop;
 
731
            betaop = 3.348 / pow(sigma[fuel[col]], 0.8189);
 
732
            A = 133 / pow(sigma[fuel[col]], 0.7913);
 
733
            gammamax =
 
734
                pow(sigma[fuel[col]],
 
735
                    1.5) / (495 + 0.0594 * pow(sigma[fuel[col]], 1.5));
 
736
            gamma =
 
737
                gammamax * pow(beta / betaop,
 
738
                               A) * exp(A * (1 - beta / betaop));
 
739
            xi = exp((0.792 + 0.681 * pow(sigma[fuel[col]], 0.5)) * (beta +
 
740
                                                                     0.1)) /
 
741
                (192 + 0.2595 * sigma[fuel[col]]);
 
742
            IR = gamma * h * (wn_dead * etaM_dead +
 
743
                              wn_live * etaM_live) * etas;
 
744
 
 
745
            R0 = IR * xi / (rhob * class_sum);
 
746
 
 
747
            if (parm.vel->answer && parm.dir->answer) {
 
748
                C = 7.47 * exp(-0.133 * pow(sigma[fuel[col]], 0.55));
 
749
                B = 0.02526 * pow(sigma[fuel[col]], 0.54);
 
750
                E = 0.715 * exp(-0.000359 * sigma[fuel[col]]);
 
751
                phiw = C * pow((double)vel[col], B) / pow(beta / betaop, E);
 
752
            }
 
753
            else
 
754
                phiw = 0.0;
 
755
 
 
756
            if (parm.slope->answer && parm.aspect->answer) {
 
757
                phis =
 
758
                    5.275 * pow(beta,
 
759
                                -0.3) * tan(0.01745 * slope[col]) *
 
760
                    tan(0.01745 * slope[col]);
 
761
            }
 
762
            else
 
763
                phis = 0.0;
 
764
 
 
765
            /*compute the maximum ROS R and its direction, 
 
766
             *vector adding for the effect of wind and slope*/
 
767
            if (parm.dir->answer && parm.aspect->answer) {
 
768
                sin_fac =
 
769
                    phiw * sin(0.01745 * dir[col]) +
 
770
                    phis * sin(0.01745 * aspect[col]);
 
771
                cos_fac =
 
772
                    phiw * cos(0.01745 * dir[col]) +
 
773
                    phis * cos(0.01745 * aspect[col]);
 
774
                phi_ws = sqrt(sin_fac * sin_fac + cos_fac * cos_fac);
 
775
                Rdir = atan2(sin_fac, cos_fac) / 0.01745;
 
776
            }
 
777
            else if (parm.dir->answer && !(parm.aspect->answer)) {
 
778
                phi_ws = phiw;
 
779
                Rdir = dir[col];
 
780
            }
 
781
            else if (!(parm.dir->answer) && parm.aspect->answer) {
 
782
                phi_ws = phis;
 
783
                Rdir = (float)aspect[col];
 
784
            }
 
785
            else {
 
786
                phi_ws = 0.0;
 
787
                Rdir = 0.0;
 
788
            }
 
789
            R = R0 * (1 + phi_ws);
 
790
            if (Rdir < 0.0)
 
791
                Rdir = Rdir + 360;
 
792
            /*printf("\nparm.dir->aanswer=%s, parm.aspect->aanswer=%s, phis=%f, phi_ws=%f, aspect[col]=%d,Rdir=%f",parm.dir->answer,parm.aspect->answer,phis,phi_ws,aspect[col],Rdir); */
 
793
 
 
794
            /*maximum spotting distance */
 
795
            if (spotting)
 
796
                spotdist[col] =
 
797
                    spot_dist(fuel[col], R, vel[col], Rdir, row, col);
 
798
 
 
799
            /*to avoid 0 R, convert ft./min to cm/min */
 
800
            R0 = 30.5 * R0;
 
801
            R = 30.5 * R;
 
802
            /*4debugging            R0 = R0/30.5/1.1; R = R/30.5/1.1; */
 
803
 
 
804
            base[col] = (int)R0;
 
805
            max[col] = (int)R;
 
806
            maxdir[col] = (int)Rdir;
 
807
            /*printf("(%d, %d)\nR0=%.2f, vel=%d, dir=%d, phiw=%.2f, s=%d, as=%d, phis=%.2f, R=%.1f, Rdir=%.0f\n", row, col, R0, vel[col], dir[col], phiw, slope[col], aspect[col], phis, R, Rdir); */
 
808
        }
 
809
        Rast_put_row(base_fd, base, CELL_TYPE);
 
810
        Rast_put_row(max_fd, max, CELL_TYPE);
 
811
        Rast_put_row(maxdir_fd, maxdir, CELL_TYPE);
 
812
        if (spotting)
 
813
            Rast_put_row(spotdist_fd, spotdist, CELL_TYPE);
 
814
    }
 
815
    G_percent(row, nrows, 2);
 
816
 
 
817
    Rast_close(fuel_fd);
 
818
    if (parm.mois_1h->answer)
 
819
        Rast_close(mois_1h_fd);
 
820
    if (parm.mois_10h->answer)
 
821
        Rast_close(mois_10h_fd);
 
822
    if (parm.mois_100h->answer)
 
823
        Rast_close(mois_100h_fd);
 
824
    Rast_close(mois_live_fd);
 
825
    if (parm.vel->answer)
 
826
        Rast_close(vel_fd);
 
827
    if (parm.dir->answer)
 
828
        Rast_close(dir_fd);
 
829
    if (parm.slope->answer)
 
830
        Rast_close(slope_fd);
 
831
    if (parm.aspect->answer)
 
832
        Rast_close(aspect_fd);
 
833
    Rast_close(base_fd);
 
834
    Rast_close(max_fd);
 
835
    Rast_close(maxdir_fd);
 
836
    if (spotting) {
 
837
        Rast_close(spotdist_fd);
 
838
        G_free(map_elev);
 
839
    }
 
840
 
 
841
    /*      
 
842
       for (model = 1; model <= 13; model++) {
 
843
       if (model == 1)
 
844
       printf("\n           Grass and Grass-dominated\n");
 
845
       else if (model == 4)
 
846
       printf("             Chaparral and Shrubfields\n");
 
847
       else if (model == 8)
 
848
       printf("                   Timber Litter\n");
 
849
       else if (model == 11)
 
850
       printf("                   Logging Slash\n");
 
851
       printf("Model %2d   ", model);
 
852
       for (class = 0; class <= 3; class++) 
 
853
       printf("%4.0f/%.3f ", SIGMA[class][model], WO[class][model]);
 
854
       printf("  %.1f  %.2f\n", DELTA[model], MX[model]);
 
855
       } 
 
856
 
 
857
       printf("\nWeight in All Fuel Subclasses:\n");
 
858
       for (model = 1; model <= 13; model++) {
 
859
       printf("model %2d  ", model);
 
860
       for (class = 0; class <= 3; class++) 
 
861
       printf("%.2f  ", Ffactor_all[class][model]);
 
862
       printf("%4.0f/%.3f=%6.0f  model %2d\n", sigma[model], wo_dead[model]+WO[3][model], sigma[model]/(wo_dead[model]+WO[3][model]), model);
 
863
       } 
 
864
       printf("\nWeight in Dead Fuel Subclasses, Dead Weight/Live Weight:\n");for (model = 1; model <= 13; model++) {
 
865
       printf("model %2d  ", model);
 
866
       for (class = 0; class <= 2; class++) 
 
867
       printf("%.2f  ", Ffactor_in_dead[class][model]);
 
868
       printf("%.2f/%.2f  model %2d\n", Ffactor_dead[model], Ffactor_live[model], model);
 
869
       } 
 
870
     */
 
871
 
 
872
    G_done_msg(_("Raster maps <%s>, <%s> and <%s> created."), name_base, name_max, name_maxdir);
 
873
    
 
874
    exit(EXIT_SUCCESS);
 
875
}