Rev 1176 | Rev 1179 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
| Rev 1176 | Rev 1178 | ||
|---|---|---|---|
| Line 361... | Line 361... | ||
| 361 | // filling of the kernel.... use the treemig
|
361 | // filling of the kernel.... use the treemig
|
| 362 | QPoint center = QPoint(kernel_offset, kernel_offset); |
362 | QPoint center = QPoint(kernel_offset, kernel_offset); |
| 363 | const float *sk_end = kernel.end(); |
363 | const float *sk_end = kernel.end(); |
| 364 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
364 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
| 365 | double d = kernel.distance(center, kernel.indexOf(p)); |
365 | double d = kernel.distance(center, kernel.indexOf(p)); |
| 366 | *p = d<=max_dist?static_cast<float>(treemig(d)):0.f; |
- | |
| - | 366 | if (d==0.) |
|
| - | 367 | *p = treemig_centercell(sqrt(cell_size*cell_size/M_PI)); // r is the radius of a circle with the same area as a cell |
|
| - | 368 | else
|
|
| - | 369 | *p = d<=max_dist?static_cast<float>(treemig(d)):0.f; |
|
| 367 | }
|
370 | }
|
| 368 | 371 | ||
| 369 | // normalize
|
372 | // normalize
|
| 370 | float sum = kernel.sum(); |
373 | float sum = kernel.sum(); |
| 371 | if (sum==0. || occupation==0.) |
374 | if (sum==0. || occupation==0.) |
| 372 | throw IException("create seed kernel: sum of probabilities = 0!"); |
375 | throw IException("create seed kernel: sum of probabilities = 0!"); |
| 373 | 376 | ||
| 374 | // the sum of all kernel cells has to equal 1
|
377 | // the sum of all kernel cells has to equal 1
|
| 375 | kernel.multiply(1.f/sum); |
- | |
| - | 378 | // kernel.multiply(1.f/sum);
|
|
| - | 379 | ||
| 376 | float fecundity_factor = static_cast<float>( max_seed / occupation); |
380 | float fecundity_factor = static_cast<float>( max_seed / occupation); |
| 377 | // probabilities are derived in multiplying by seed number, and dividing by occupancy criterion
|
381 | // probabilities are derived in multiplying by seed number, and dividing by occupancy criterion
|
| 378 | kernel.multiply( fecundity_factor ); |
382 | kernel.multiply( fecundity_factor ); |
| 379 | // all cells that get more seeds than the occupancy criterion are considered to have no seed limitation for regeneration
|
383 | // all cells that get more seeds than the occupancy criterion are considered to have no seed limitation for regeneration
|
| 380 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
384 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
| 381 | *p = qMin(*p, 1.f); |
385 | *p = qMin(*p, 1.f); |
| 382 | }
|
386 | }
|
| 383 | // set the parent cell to 1
|
387 | // set the parent cell to 1
|
| 384 | kernel.valueAtIndex(kernel_offset, kernel_offset)=1.f; |
- | |
| - | 388 | //kernel.valueAtIndex(kernel_offset, kernel_offset)=1.f;
|
|
| 385 | 389 | ||
| 386 | 390 | ||
| 387 | // some final statistics....
|
391 | // some final statistics....
|
| 388 | if (logLevelInfo()) qDebug() << "kernel setup.Species:"<< mSpecies->id() << "kernel-size: " << kernel.sizeX() << "x" << kernel.sizeY() << "pixels, sum: " << kernel.sum(); |
392 | if (logLevelInfo()) qDebug() << "kernel setup.Species:"<< mSpecies->id() << "kernel-size: " << kernel.sizeX() << "x" << kernel.sizeY() << "pixels, sum: " << kernel.sum(); |
| 389 | 393 | ||
| Line 408... | Line 412... | ||
| 408 | mLDDDistance.push_back(mLDDDistance.last() + (r_max-r_min)/static_cast<float>(mLDDRings)); |
412 | mLDDDistance.push_back(mLDDDistance.last() + (r_max-r_min)/static_cast<float>(mLDDRings)); |
| 409 | double r_out = mLDDDistance.last(); |
413 | double r_out = mLDDDistance.last(); |
| 410 | // calculate the value of the kernel for the middle of the ring
|
414 | // calculate the value of the kernel for the middle of the ring
|
| 411 | double ring_in = treemig(r_in); // kernel value at the inner border of the ring |
415 | double ring_in = treemig(r_in); // kernel value at the inner border of the ring |
| 412 | double ring_out = treemig(r_out); // kernel value at the outer border of the ring |
416 | double ring_out = treemig(r_out); // kernel value at the outer border of the ring |
| 413 | double ring_val = ring_in - ring_out; // this p should be distributed within the ring |
- | |
| - | 417 | double ring_val = (ring_in + ring_out)/2.; // this is the average p |
|
| 414 | ring_val *= mFecundityFactor; // include fecundity (along the lines of the kernel calculations) |
418 | ring_val *= mFecundityFactor; // include fecundity (along the lines of the kernel calculations) |
| 415 | // calculate the area of the ring
|
419 | // calculate the area of the ring
|
| 416 | double ring_area = (r_out*r_out - r_in*r_in)*M_PI; // in square meters |
420 | double ring_area = (r_out*r_out - r_in*r_in)*M_PI; // in square meters |
| 417 | double ring_px = ring_area / (mSeedMap.cellsize()*mSeedMap.cellsize()); // # of seed pixels on the ring |
421 | double ring_px = ring_area / (mSeedMap.cellsize()*mSeedMap.cellsize()); // # of seed pixels on the ring |
| 418 | double n_px = ring_val * ring_px / mLDDProbability; |
422 | double n_px = ring_val * ring_px / mLDDProbability; |
| Line 433... | Line 437... | ||
| 433 | }
|
437 | }
|
| 434 | */
|
438 | */
|
| 435 | 439 | ||
| 436 | /// the used kernel function
|
440 | /// the used kernel function
|
| 437 | /// see also Appendix B of iland paper II (note the different variable names)
|
441 | /// see also Appendix B of iland paper II (note the different variable names)
|
| - | 442 | /// the function returns the seed density at a point with distance 'distance'.
|
|
| 438 | double SeedDispersal::treemig(const double &distance) |
443 | double SeedDispersal::treemig(const double &distance) |
| 439 | {
|
444 | {
|
| 440 | double p1 = (1.-mTM_ks)*exp(-distance/mTM_as1)/mTM_as1; |
445 | double p1 = (1.-mTM_ks)*exp(-distance/mTM_as1)/mTM_as1; |
| 441 | double p2 = 0.; |
446 | double p2 = 0.; |
| 442 | if (mTM_as2>0.) |
447 | if (mTM_as2>0.) |
| 443 | p2 = mTM_ks*exp(-distance/mTM_as2)/mTM_as2; |
448 | p2 = mTM_ks*exp(-distance/mTM_as2)/mTM_as2; |
| 444 | return p1 + p2; |
- | |
| - | 449 | double s = p1 + p2; |
|
| - | 450 | // 's' is the density for radius 'distance' - not for specific point with that distance.
|
|
| - | 451 | // (i.e. the integral over the one-dimensional treemig function is 1, but if applied for 2d cells, the
|
|
| - | 452 | // sum would be much larger as all seeds arriving at 'distance' would be arriving somewhere at the circle with radius 'distance')
|
|
| - | 453 | // convert that to a density at a point, by dividing with the circumference at the circle with radius 'distance'
|
|
| - | 454 | s = s / (2.*std::max(distance, 0.01)*M_PI); |
|
| - | 455 | ||
| - | 456 | return s; |
|
| - | 457 | }
|
|
| - | 458 | ||
| - | 459 | double SeedDispersal::treemig_centercell(const double &max_distance) |
|
| - | 460 | {
|
|
| - | 461 | // use 100 steps and calculate dispersal kernel for consecutive rings
|
|
| - | 462 | double sum = 0.; |
|
| - | 463 | for (int i=0;i<100;i++) { |
|
| - | 464 | double r_in = i*max_distance/100.; |
|
| - | 465 | double r_out = (i+1)*max_distance/100.; |
|
| - | 466 | double ring_area = (r_out*r_out-r_in*r_in)*M_PI; |
|
| - | 467 | // the value of each ring is: treemig(r) * area of the ring
|
|
| - | 468 | sum += treemig((r_out+r_in)/2.)*ring_area; |
|
| - | 469 | }
|
|
| - | 470 | return sum; |
|
| 445 | }
|
471 | }
|
| 446 | 472 | ||
| 447 | /// calculate the distance where the probability falls below 'value'
|
473 | /// calculate the distance where the probability falls below 'value'
|
| 448 | double SeedDispersal::treemig_distanceTo(const double value) |
474 | double SeedDispersal::treemig_distanceTo(const double value) |
| 449 | {
|
475 | {
|