WebSVN - public iLand - Diff - Rev 1181 and 1182

+            }
             *p=1.f; // mark as processed
         } // *p==1
     } // for()
+}
+// because C modulo operation gives negative numbers for negative values, here a fix
+// that always returns positive numbers: http://www.lemoda.net/c/modulo-operator/
+#define MOD(a,b) ((((a)%(b))+(b))%(b))
 void SeedDispersal::distributeSeeds(Grid<float> *seed_map)
+{
     Grid<float> &sourcemap = seed_map ? *seed_map : mSourceMap; // switch to extra seed map if provided
     Grid<float> &kernel = mKernelSeedYear;
+    // *** estimate seed production (based on leaf area) ***
     // calculate number of seeds; the source map holds now m2 leaf area on 20x20m pixels
     // after this step, each source cell has a value between 0 (no source) and 1 (fully covered cell)
     float fec = species()->fecundity_m2();
     if (!species()->isSeedYear())
         fec *= mNonSeedYearFraction;
-//    mSeedMap.initialize(0.f); // just for debugging...
+    //    mSeedMap.initialize(0.f); // just for debugging...
     int offset = kernel.sizeX() / 2; // offset is the index of the center pixel
     // source mode
+    // *** seed distribution (Kernel + long distance dispersal) ***
+    if (GlobalSettings::instance()->model()->settings().torusMode==false) {
+        // ** standard case (no torus) **
-    for (float *src=sourcemap.begin(); src!=sourcemap.end(); ++src) {
+        for (float *src=sourcemap.begin(); src!=sourcemap.end(); ++src) {
-        if (*src>0.f) {
+            if (*src>0.f) {
-            QPoint sm=sourcemap.indexOf(src)-QPoint(offset, offset);
+                QPoint sm=sourcemap.indexOf(src)-QPoint(offset, offset);
-            int sx = sm.x(), sy=sm.y();
+                int sx = sm.x(), sy=sm.y();
-            for (int iy=0;iy<kernel.sizeY();++iy) {
+                for (int iy=0;iy<kernel.sizeY();++iy) {
+                    }
+                }
+            }
+        }
+    } else {
+        // **** seed distribution in torus mode ***
+        int seedmap_offset = sourcemap.indexAt(QPointF(0., 0.)).x(); // the seed maps have x extra rows/columns
+        QPoint torus_pos;
+        int seedpx_per_ru = static_cast<int>((cRUSize/sourcemap.cellsize()));
+        for (float *src=sourcemap.begin(); src!=sourcemap.end(); ++src) {
+            if (*src>0.f) {
+                QPoint sm=sourcemap.indexOf(src);
+                // get the origin of the resource unit *on* the seedmap in *seedmap-coords*:
+                QPoint offset_ru( ((sm.x()-seedmap_offset) / seedpx_per_ru) * seedpx_per_ru + seedmap_offset,
+                                 ((sm.y()-seedmap_offset) / seedpx_per_ru) * seedpx_per_ru + seedmap_offset);  // coords RU origin
+                QPoint offset_in_ru((sm.x()-seedmap_offset) % seedpx_per_ru, (sm.y()-seedmap_offset) % seedpx_per_ru );  // offset of current point within the RU
+                //QPoint sm=sourcemap.indexOf(src)-QPoint(offset, offset);
+                for (int iy=0;iy<kernel.sizeY();++iy) {
+                    for (int ix=0;ix<kernel.sizeX();++ix) {
+                        torus_pos = offset_ru + QPoint(MOD((offset_in_ru.x() - offset + ix), seedpx_per_ru), MOD((offset_in_ru.y() - offset + iy), seedpx_per_ru));
+                        if (mSeedMap.isIndexValid(torus_pos))
+                            mSeedMap.valueAtIndex(torus_pos)+= *src * kernel(ix, iy);
+                    }
+                }
+                // long distance dispersal
+                if (!mLDDDensity.isEmpty()) {
+                    for (int r=0;r<mLDDDensity.size(); ++r) {
+                        float ldd_val = mLDDSeedlings / fec; // pixels will have this probability [note: fecundity will be multiplied below]
+                        int n;
+                        if (mLDDDensity[r]<1)
+                            n = drandom()<mLDDDensity[r] ? 1 : 0;
+                        else
+                            n = round( mLDDDensity[r] ); // number of pixels to activate
+                        for (int i=0;i<n;++i) {
+                            // distance and direction:
+                            double radius = nrandom(mLDDDistance[r], mLDDDistance[r+1]) / mSeedMap.cellsize(); // choose a random distance (in pixels)
+                            double phi = drandom()*2.*M_PI; // choose a random direction
+                            QPoint ldd( radius*cos(phi),  + radius*sin(phi)); // destination (offset)
+                            torus_pos = offset_ru + QPoint(MOD((offset_in_ru.x()+ldd.x()),seedpx_per_ru), MOD((offset_in_ru.y()+ldd.y()),seedpx_per_ru) );
+                            if (mSeedMap.isIndexValid(torus_pos)) {
+                                float &val = mSeedMap.valueAtIndex(torus_pos);
+                                _debug_ldd++;
+                                val += ldd_val;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    } // torus
     // now the seed sources (0..1) are spatially distributed by the kernel (and LDD) without altering the magnitude;
     // now we include the fecundity (=seedling potential per m2 crown area), and convert to the establishment probability p_seed.

Subversion Repositories public iLand

(root)/src/core/seeddispersal.cpp @ 1222 - Rev 1181 → 1182