Logo Search packages:      
Sourcecode: saga version File versions  Download package

Flow_Parallel.cpp


///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                    Module Library:                    //
//                     ta_hydrology                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                   Flow_Parallel.cpp                   //
//                                                       //
//                 Copyright (C) 2003 by                 //
//                      Olaf Conrad                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// This file is part of 'SAGA - System for Automated     //
// Geoscientific Analyses'. SAGA is free software; you   //
// can redistribute it and/or modify it under the terms  //
// of the GNU General Public License as published by the //
// Free Software Foundation; version 2 of the License.   //
//                                                       //
// SAGA is distributed in the hope that it will be       //
// useful, but WITHOUT ANY WARRANTY; without even the    //
// implied warranty of MERCHANTABILITY or FITNESS FOR A  //
// PARTICULAR PURPOSE. See the GNU General Public        //
// License for more details.                             //
//                                                       //
// You should have received a copy of the GNU General    //
// Public License along with this program; if not,       //
// write to the Free Software Foundation, Inc.,          //
// 59 Temple Place - Suite 330, Boston, MA 02111-1307,   //
// USA.                                                  //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     oconrad@saga-gis.org                   //
//                                                       //
//    contact:    Olaf Conrad                            //
//                Institute of Geography                 //
//                University of Goettingen               //
//                Goldschmidtstr. 5                      //
//                37077 Goettingen                       //
//                Germany                                //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "Flow_Parallel.h"


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CFlow_Parallel::CFlow_Parallel(void)
{
      Set_Name          (_TL("Parallel Processing"));

      Set_Author        (SG_T("(c) 2001 by O.Conrad"));

      Set_Description   (_TW(
            "Parallel processing of cells for calculation of flow accumulation and related parameters. "
            "This set of algorithms processes a DEM downwards from the highest to the lowest cell.\n\n"

            "References:\n\n"

            "Deterministic 8\n"
            "- O'Callaghan, J.F. / Mark, D.M. (1984):\n"
            "    'The extraction of drainage networks from digital elevation data',\n"
            "    Computer Vision, Graphics and Image Processing, 28:323-344\n\n"

            "Rho 8:\n"
            "- Fairfield, J. / Leymarie, P. (1991):\n"
            "    'Drainage networks from grid digital elevation models',\n"
            "    Water Resources Research, 27:709-717\n\n"

            "Braunschweiger Reliefmodell:\n"
            "- Bauer, J. / Rohdenburg, H. / Bork, H.-R. (1985):\n"
            "    'Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse',\n"
            "    Landschaftsgenese und Landschaftsoekologie, H.10, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle,\n"
            "    Grundlagenarbeiten zu Analyse von Agrar-Oekosystemen, (Eds.: Bork, H.-R. / Rohdenburg, H.), p.1-15\n\n"

            "Deterministic Infinity:\n"
            "- Tarboton, D.G. (1997):\n"
            "    'A new method for the determination of flow directions and upslope areas in grid digital elevation models',\n"
            "    Water Ressources Research, Vol.33, No.2, p.309-319\n\n"

            "Multiple Flow Direction:\n"
            "- Freeman, G.T. (1991):\n"
            "    'Calculating catchment area with divergent flow based on a regular grid',\n"
            "    Computers and Geosciences, 17:413-22\n\n"

            "- Quinn, P.F. / Beven, K.J. / Chevallier, P. / Planchon, O. (1991):\n"
            "    'The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models',\n"
            "    Hydrological Processes, 5:59-79\n\n"
            
            "Triangular Multiple Flow Direction\n"
            "- Seibert, J. / McGlynn, B. (2007):\n"
            "    'A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models',\n"
            "    Water Ressources Research, Vol. 43, W04501\n"
            "    C++ Implementation into SAGA by Thomas Grabs, Copyrights (c) 2007\n"
            "    Contact: thomas.grabs@natgeo.su.se, jan.seibert@natgeo.su.se \n\n"
      ));


      //-----------------------------------------------------
      // Output...

      Parameters.Add_Grid(
            NULL  , "CASPECT"       , _TL("Catchment Aspect"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "FLWPATH"       , _TL("Flow Path Length"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );


      //-----------------------------------------------------
      // Method...

      Parameters.Add_Choice(
            NULL  , "Method"        , _TL("Method"),
            _TL(""),
            CSG_String::Format(SG_T("%s|%s|%s|%s|%s|%s|"),
                  _TL("Deterministic 8"),
                  _TL("Rho 8"),
                  _TL("Braunschweiger Reliefmodell"),
                  _TL("Deterministic Infinity"),
                  _TL("Multiple Flow Direction"),
                  _TL("Multiple Triangular Flow Directon")
            ), 4
      );


      //-----------------------------------------------------
      // Options...

      Parameters.Add_Value(
            NULL  , "DOLINEAR"      , _TL("Linear Flow"),
            _TL("Use D8 if catchment area becomes higher than specified threshold."),
            PARAMETER_TYPE_Bool
      );

      Parameters.Add_Value(
            NULL  , "LINEARTHRS"    , _TL("Linear Flow Threshold"),
            _TL("Use D8 if catchment area becomes higher than specified threshold (Cells)."),
            PARAMETER_TYPE_Double   ,     500
      );

      Parameters.Add_Value(
            NULL  , "CONVERGENCE"   , _TL("Convergence"),
            _TL("Convergence factor for Multiple Flow Direction Algorithm after Freeman (1991) and Multiple Triangular Flow Directon Algorithm after Seibert/Glynn (2007)."),
            PARAMETER_TYPE_Double   , 1.1, 0.0, true
      );
}

//---------------------------------------------------------
CFlow_Parallel::~CFlow_Parallel(void)
{}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::On_Initialize(void)
{
      pCatch_Aspect     = Parameters("CASPECT")->asGrid();
      pFlowPath         = Parameters("FLWPATH")->asGrid();

      TH_LinearFlow     = Parameters("DOLINEAR")->asBool() && pDTM
                              ? Parameters("LINEARTHRS")->asDouble() // (pDTM->Get_NX() * pDTM->Get_NY())
                              : -1.0;

      MFD_Converge      = Parameters("CONVERGENCE")->asDouble();
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Calculate(void)
{
      int         x, y;

      for(y=0; y<Get_NY() && Set_Progress(y); y+=Step)
      {
            for(x=0; x<Get_NX(); x+=Step)
            {
                  Init_Cell(x, y);
            }
      }

      Set_Flow();
}

//---------------------------------------------------------
void CFlow_Parallel::Calculate(int x, int y)
{
      Init_Cell(x, y);

      Set_Flow();
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_Flow(void)
{
      int         Method, x, y;

      long  n;

      //-----------------------------------------------------
      Method      = Parameters("Method")->asInt();

      if( Method == 2 )
      {
            BRM_Init();
      }

      //-----------------------------------------------------
      for(n=0; n<Get_NCells() && Set_Progress_NCells(n); n++)
      {
            pDTM->Get_Sorted(n,x,y, true, false);

            //if( !(n%(4*Get_NX())) )DataObject_Update(pFlow);

            if( TH_LinearFlow > 0.0 && pCatch->asDouble(x, y) >= TH_LinearFlow )
            {
                  Set_D8(x, y);
            }
            else
            {
                  switch( Method )
                  {
                  case 0:           Set_D8            (x, y);           break;
                  case 1:           Set_Rho8    (x, y);           break;
                  case 3:           Set_DInf    (x, y);           break;
                  case 4:           Set_MFD           (x, y);           break;
                  case 5:           Set_MDInf   (x, y);           break;                        
                  case 2:           Set_BRM           (x, y);           break;
                  }
            }
      }

      //-----------------------------------------------------
      if( pRoute )
      {
            for(n=0; n<Get_NCells() && Set_Progress_NCells(n); n++)
            {
                  if( pDTM->Get_Sorted(n, x, y, false) )
                  {
                        Check_Route(x, y);
                  }
            }
      }
}

//---------------------------------------------------------
void CFlow_Parallel::Check_Route(int x, int y)
{
      bool  bSink;
      int         i, ix, iy;
      double      z;

      if( pRoute->asChar(x, y) > 0 )
      {
            z     = pDTM->asDouble(x, y);

            for(i=0, bSink=true; i<8 && bSink; i++)
            {
                  ix    = Get_xTo(i, x);
                  iy    = Get_yTo(i, y);

                  if( !is_InGrid(ix, iy) || z > pDTM->asDouble(ix, iy) )
                  {
                        bSink = false;
                  }
            }

            //-------------------------------------------------
            if( bSink )
            {
                  i     = pRoute->asChar(x, y);

                  ix    = Get_xTo(i, ix);
                  iy    = Get_yTo(i, iy);

                  //---------------------------------------------
                  while( is_InGrid(ix, iy) )
                  {
                        Add_Portion(x, y, ix, iy);

                        i     = pRoute->asChar(ix, iy);

                        if( i > 0 )
                        {
                              ix    = Get_xTo(i, ix);
                              iy    = Get_yTo(i, iy);
                        }
                        else
                        {
                              i     = pDTM->Get_Gradient_NeighborDir(ix, iy);

                              if( i >= 0 )
                              {
                                    ix    = Get_xTo(i, ix);
                                    iy    = Get_yTo(i, iy);
                              }
                              else
                              {
                                    ix    = -1;
                              }
                        }
                  }
            }
      }
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                            Deterministic 8                                  //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_D8(        int x, int y )
{
      int         Direction;

      Direction   = pDTM->Get_Gradient_NeighborDir(x, y);

      if( Direction >= 0 )
      {
            Add_Fraction(x, y, Direction);
      }
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                  Rho 8                                      //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_Rho8(      int x, int y )
{
      int         i, ix, iy, iMax;

      double      z, d, dMax;

      z           = pDTM->asDouble(x, y);
      iMax  = -1;

      for(i=0; i<8; i++)
      {
            ix          = Get_xTo(i, x);
            iy          = Get_yTo(i, y);

            if( !pDTM->is_InGrid(ix, iy) )
            {
                  return;
            }
            else
            {
                  d           = z - pDTM->asDouble(ix, iy);

                  if( i % 2 == 1 )
                  {
                        d           /= 1.0 + rand() / (double)RAND_MAX;
                  }

                  if( iMax < 0 || (iMax >= 0 && d > dMax) )
                  {
                        iMax  = i;
                        dMax  = d;
                  }
            }
      }

      if( iMax >= 0 )
      {
            Add_Fraction(x, y, iMax);
      }
}

/*void CFlow_Parallel::Set_Rho8(    int x, int y )
{
      int         Direction, ix, iy;

      double      Slope, Aspect;

      Get_Gradient(x, y, Slope, Aspect);

      if( Aspect >= 0 )
      {
            Direction   = (int)(Aspect / M_PI_045);

            if( fmod(Aspect, M_PI_045) / M_PI_045 > rand() / (double)RAND_MAX )
            {
                  Direction++;
            }

            Direction   %= 8;

            ix                = Get_xTo(Direction, x);
            iy                = Get_yTo(Direction, y);

            if( is_InGrid(ix, iy) && pDTM->asDouble(ix, iy) >= pDTM->asDouble(x, y) )
            {
                  Direction   = Get_Direction_Lowest(pDTM, x, y);
            }

            Add_Fraction(x, y, Direction);
      }
}/**/


///////////////////////////////////////////////////////////
//                                                                                   //
//                      Deterministic Infinity                           //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_DInf(      int x, int y )
{
      int         Direction, x1, y1, x2, y2;

      double      Slope, Aspect, z;

      Get_Gradient(x, y, Slope, Aspect);

      if( Aspect >= 0.0 )
      {
            Direction   = (int)(Aspect / M_PI_045);
            Aspect            =  fmod(Aspect , M_PI_045) / M_PI_045;

            z                 = pDTM->asDouble(x, y);
            x1                = Get_xTo(Direction + 0, x);
            y1                = Get_yTo(Direction + 0, y);
            x2                = Get_xTo(Direction + 1, x);
            y2                = Get_yTo(Direction + 1, y);

            if( (!is_InGrid(x1, y1) || z > pDTM->asDouble(x1, y1))
            &&    (!is_InGrid(x2, y2) || z > pDTM->asDouble(x2, y2)) )
            {
                  Add_Fraction(x, y,  Direction             , 1.0 - Aspect);
                  Add_Fraction(x, y, (Direction + 1) % 8    ,       Aspect);
            }
            else
            {
                  Add_Fraction(x, y, pDTM->Get_Gradient_NeighborDir(x, y));
            }
      }
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                      Multiple Flow Direction                          //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_MFD( int x, int y )
{
      int         i, ix, iy;

      double      z, d, dzSum, dz[8];

      z           = pDTM->asDouble(x, y);
      dzSum = 0.0;

      for(i=0; i<8; i++)
      {
            ix          = Get_xTo(i, x);
            iy          = Get_yTo(i, y);

            if( pDTM->is_InGrid(ix, iy) )
            {
                  d           = z - pDTM->asDouble(ix, iy);
            }
            else
            {
                  ix          = Get_xTo(i + 4, x);
                  iy          = Get_yTo(i + 4, y);

                  if( pDTM->is_InGrid(ix, iy) )
                  {
                        d           = pDTM->asDouble(ix, iy) - z;
                  }
                  else
                  {
                        d           = 0.0;
                  }
            }

            if( d > 0.0 )
            {
                  dzSum += (dz[i]   = pow(d / Get_Length(i), MFD_Converge));
            }
            else
            {
                  dz[i] = 0.0;
            }
      }

      if( dzSum )
      {
            for(i=0; i<8; i++)
            {
                  if( dz[i] )
                  {
                        Add_Fraction(x, y, i, dz[i] / dzSum );
                  }
            }
      }
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                Multiple Triangular Flow Directon                //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_MDInf(     int x, int y )
{
      int         i, ii, ix, iy;
      double      z, dzSum, dz[8];
      double      hs, hr, s_facet[8], r_facet[8];
      double      valley[8], portion[8];
      double  cellsize, cellarea;
      double      nx, ny, nz, n_norm;
      bool  Dir_inGrid[8];

      z                 = pDTM->asDouble(x, y);
      cellsize    = Get_Cellsize();
      cellarea    = cellsize * cellsize;

      for(i=0; i<8; i++)
      {
            ix          = Get_xTo(i, x);
            iy          = Get_yTo(i, y);
            
//    ???
//          if( x == 6 && y == 7 )
//          {
//                int k = 1;
//          }
//    ???

            s_facet[i] = r_facet[i] = -999.0;

            if( pDTM->is_InGrid(ix, iy) )
            {
                  dz[i]             = z - pDTM->asDouble(ix, iy);
                  Dir_inGrid[i]     = true;
            }
            else
            {
                  dz[i]             = 0.0;
                  Dir_inGrid[i]     = false;
            }
      }

      for(i=0; i<8; i++)
      {
            hr = hs = -999.0;

            if( Dir_inGrid[i] )
            {
                  ii = (i < 7) ? i+1 : 0;

                  if( Dir_inGrid[ii] )
                  {
                        nx = ( dz[ii] * Get_yTo( i  ) - dz[i ] * Get_yTo( ii )) * cellsize;
                        ny = ( dz[i ] * Get_xTo( ii ) - dz[ii] * Get_xTo( i  )) * cellsize;

                        nz = ( Get_xTo( i ) * Get_yTo( ii ) - Get_xTo( ii ) * Get_yTo( i )) * cellarea;

                        n_norm = sqrt( nx*nx + ny*ny +nz*nz );

                        if( nx == 0.0 )
                        {
                              hr = (ny >= 0.0)? 0.0 : M_PI;
                        } 
                        else if( nx < 0.0 )
                        {
                              hr = M_PI_270 - atan(ny / nx);
                        } 
                        else
                        {
                              hr = M_PI_090 - atan(ny / nx);
                        }

                        hs = -tan( acos( nz/n_norm ) );

                        if( hr < i * M_PI_045 || hr > (i+1) * M_PI_045 )
                        {
                              if( dz[i] > dz[ii] )
                              {
                                    hr = i * M_PI_045;
                                    hs = dz[i] / Get_Length(i);
                              }
                              else
                              {
                                    hr = ii * M_PI_045;
                                    hs = dz[ii] / Get_Length(ii);                               
                              }
                        }
                        
                  }
                  else if( dz[i] > 0.0 )
                  {
                        hr = i * M_PI_045;
                        hs = dz[i] / Get_Length(i);
                  }
                  
                  s_facet[i] = hs;
                  r_facet[i] = hr;
            }
      }
      
      dzSum       = 0.0;
      
      for(i=0; i<8; i++)
      {           
            valley[i]   = 0.0;
            ii = (i < 7)? i+1 : 0;
            
            if( s_facet[i] > 0.0 )
            {
                  if( r_facet[i] > i * M_PI_045 && r_facet[i] < (i+1) * M_PI_045 )
                  {
                        valley[i] = s_facet[i];
                  }
                  else if( r_facet[i] == r_facet[ii] )
                  {
                        valley[i] = s_facet[i];
                  }
                  else if( s_facet[ii] == -999.0 && r_facet[i] == (i+1) * M_PI_045)
                  {
                        valley[i] = s_facet[i];
                  }
                  else
                  {
                        ii = (i > 0)? i-1 : 7;

                        if( s_facet[ii] == -999.0 && r_facet[i] == i * M_PI_045 )
                        {
                              valley[i] = s_facet[i];
                        }
                  }
                  
                  valley[i] = pow(valley[i], MFD_Converge);
                  dzSum += valley[i];
            } 
            
            portion[i] = 0.0;
      }

      if( dzSum )
      {
            for(i=0; i<8; i++)
            {
                  if (i < 7)
                  {
                        ii = i+1;
                  }
                  else
                  {
                        ii = 0;

                        if( r_facet[i] == 0.0)
                        {
                              r_facet[i] = M_PI_360;
                        }
                  }

                  if( valley[i] )
                  {
                        valley[i] /= dzSum;

                        portion[i] += valley[i] * ((i+1) * M_PI_045 - r_facet[i]) / M_PI_045;

                        portion[ii] += valley[i] * (r_facet[i] - i * M_PI_045) / M_PI_045;
                  }
            }

            for(i=0; i<8; i++)
            {
                  Add_Fraction(x, y, i, portion[i]);
            }
      }
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                      Braunschweiger Reliefmodell                      //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "Flow_BRM.h"

#define BRM_nint(x)           ((int)(x >= 0 ? x + 0.5 : x - 0.5))


//---------------------------------------------------------
void CFlow_Parallel::Set_BRM( int x, int y )
{
      int         Dir, QBinaer,
                  ix[3], iy[3],
                  nexp[6];

      double      QLinks, QMitte, QRecht,
                  nnei[6];

      //-----------------------------------------------------
      if( is_InGrid(x, y, 1) ) // Rand !!!
      {
            Dir   = BRM_InitRZ(x,y,ix,iy);

            if( Dir >= 0 )
            {
                  if( Dir % 2 )
                  {
                        BRM_GetDiago(Dir,x,y,ix,iy,nnei,nexp);
                        BRM_QStreuung(4,1,nnei,nexp,QBinaer,QLinks,QMitte,QRecht);
                  }
                  else
                  {
                        BRM_GetOrtho(Dir,x,y,ix,iy,nnei,nexp);
                        BRM_QStreuung(6,0,nnei,nexp,QBinaer,QLinks,QMitte,QRecht);
                  }

                  Add_Fraction(x,y,(Dir+1)%8,BRM_BitMtrx[0][QBinaer] ? QLinks : 0);
                  Add_Fraction(x,y,(Dir+0)%8,BRM_BitMtrx[1][QBinaer] ? QMitte : 0);
                  Add_Fraction(x,y,(Dir+7)%8,BRM_BitMtrx[2][QBinaer] ? QRecht : 0);
            }
      }
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_Init(void)
{
      int   i;

      double      DXT   = Get_Cellsize()/2,
                  DYT   = Get_Cellsize()/2;

      //-----------------------------------------------------
      BRM_kgexp[0]      = (int)(atan2(DXT     , Get_Cellsize()) * M_RAD_TO_DEG);
      BRM_kgexp[1]      = (int)(atan2(Get_Cellsize(), DYT     ) * M_RAD_TO_DEG) + 1;
      BRM_kgexp[2]      = (int)(atan2(Get_Cellsize(),-DYT     ) * M_RAD_TO_DEG);
      BRM_kgexp[3]      = (int)(atan2(DXT     ,-Get_Cellsize()) * M_RAD_TO_DEG) + 1;

      for(i=0; i<4; i++)
            BRM_kgexp[i+4]    = BRM_kgexp[i] + 180;

      //-----------------------------------------------------
      for(i=0; i<=360; i++)
      {
            BRM_sinus[i]      = -sin(i * M_DEG_TO_RAD);
            BRM_cosin[i]      = -cos(i * M_DEG_TO_RAD);
      }

      //---BRM_idreh---------------------------------------------
      BRM_idreh[0]      = 180;
      BRM_idreh[1]      = 180 -     BRM_nint(atan2(Get_Cellsize(), Get_Cellsize()) * M_RAD_TO_DEG);
      BRM_idreh[2]      =  90;
      BRM_idreh[3]      =           BRM_nint(atan2(Get_Cellsize(), Get_Cellsize()) * M_RAD_TO_DEG);
      BRM_idreh[4]      =   0;

      for(i=1; i<4; i++)
            BRM_idreh[i+4]    = BRM_idreh[i] + 180;
}

//---------------------------------------------------------
int CFlow_Parallel::BRM_InitRZ(int x, int y, int ix[3], int iy[3])
{
      int         i, j, Dir;

      double      Slope, Aspect;

      Get_Gradient(x, y, Slope, Aspect);

      Aspect      *= M_RAD_TO_DEG;

      if( Aspect < 0 )
      {
            return( -1 );
      }

      //---Kategorisierte-Exposition-------------------------
      Dir   = 0;
      while( Aspect > BRM_kgexp[Dir] && Dir < 8 )
            Dir++;
      Dir   %=8;

    //---Finde-Die-3-ZielRasterZellen----------------------
      for(i=0; i<3; i++)      // zxy[]: 0=Recht, 1=Mitte, 2=Links
      {
            j           = (Dir + 7 + i) % 8;
            ix[2-i]     = Get_xTo(j,x);
            iy[2-i]     = Get_yTo(j,y);
      }

      return(Dir);
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_GetOrtho(int Dir, int x, int y, int ix[3], int iy[3], double nnei[6], int nexp[6])
{
      int         jx, jy, i,
                  i0    = (Dir + 2) % 8,
                  i4    = (Dir + 6) % 8;

      double      Slope, Aspect;

      for(i=0; i<3; i++)
      {
            jx          = ix[i];
            jy          = iy[i];

            Get_Gradient(jx, jy, Slope, Aspect);

            nnei[i]     = M_RAD_TO_DEG * Slope;
            nexp[i]     = (int)(M_RAD_TO_DEG * Aspect);
      }

      jx          = Get_xTo(i0,x);
      jy          = Get_yTo(i0,y);

      Get_Gradient(jx, jy, Slope, Aspect);

      nnei[3]     = M_RAD_TO_DEG * Slope;
      nexp[3]     = (int)(M_RAD_TO_DEG * Aspect);

      jx          = Get_xTo(i4,x);
      jy          = Get_yTo(i4,y);

      Get_Gradient(jx, jy, Slope, Aspect);

      nnei[5]     = M_RAD_TO_DEG * Slope;
      nexp[5]     = (int)(M_RAD_TO_DEG * Aspect);

      Get_Gradient(x, y, Slope, Aspect);

      nnei[4]     = M_RAD_TO_DEG * Slope;
      nexp[4]     = (int)(M_RAD_TO_DEG * Aspect);     //[jy][jx]) ????!!!!...;

      for(i=0; i<6; i++)
            if(nexp[i]<0)
                  nexp[i]     = nexp[4];

      for(i=0; i<6; i++)
      {
            nexp[i]     += BRM_idreh[Dir];

            if(nexp[i]>360)
                  nexp[i]     -= 360;
      }
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_GetDiago(int Dir, int x, int y, int ix[3], int iy[3], double nnei[6], int nexp[6])
{
      int         i;

      double      Slope, Aspect;

      Get_Gradient( x   ,  y   , Slope, Aspect);
      nexp[0]     = (int)(M_RAD_TO_DEG * Aspect);
      nnei[0]     = M_RAD_TO_DEG * Slope;

      Get_Gradient(ix[0], iy[0], Slope, Aspect);
      nexp[1]     = (int)(M_RAD_TO_DEG * Aspect);
      nnei[1]     = M_RAD_TO_DEG * Slope;

      Get_Gradient(ix[2], iy[2], Slope, Aspect);
      nexp[2]     = (int)(M_RAD_TO_DEG * Aspect);
      nnei[2]     = M_RAD_TO_DEG * Slope;

      Get_Gradient(ix[1], iy[1], Slope, Aspect);
      nexp[3]     = (int)(M_RAD_TO_DEG * Aspect);
      nnei[3]     = M_RAD_TO_DEG * Slope;

      for(i=1; i<4; i++)
            if(nexp[i]<0)
                  nexp[i]     = nexp[0];

      for(i=0; i<4; i++)
      {
            nexp[i]     += BRM_idreh[Dir];

            if(nexp[i]>360)
                  nexp[i]     -= 360;
      }
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_QStreuung(int i64, int g64, double nnei[6], int nexp[6], int &QBinaer, double &QLinks, double &QMitte, double &QRecht)
{
      int         i, j, ix, iy,
                  ALinks, AMitte=2, ARecht;

      double      x=1, y=1, sg=0, a,
                  s[6], c[6];

      ALinks      = ARecht    = 0;
      QLinks      = QRecht    = 0.0;

      for(i=0; i<i64; i++)
            sg    += nnei[i];

      sg    = i64 / sg;

      for(i=0; i<i64; i++)
      {
            a           = sg * nnei[i];
            j           = nexp[i];
            s[i]  = a * BRM_sinus[j];
            c[i]  = a * BRM_cosin[j];
      }

    //---QLinks-ermitteln----------------------------------
      for(i=0; i<100; i++)
      {
            ix    = BRM_nint(x) - 1;
            iy    = BRM_nint(y) - 1;

            for(j=0; j<i64; j++)
            {
                  a           = BRM_g[g64][j][ix][iy];
                  x           += s[j] * a;
                  y           += c[j] * a;
            }      

            if(x<1)
            {
                  ALinks      = 0;
                  QLinks      = 0;
                  break;
            }

            if( x>8.99 || y<1 )
            {
                  ALinks      = 4;
                  QLinks      = 1;              
                  break;
            }

            if(y>8.95)
            {
                  if(x<1.02)
                  {
                        ALinks      = 0;
                        QLinks      = 0;
                  }
                  else
                  {
                        ALinks      = 4;
                        QLinks      = (x - 1) / 8;

                        if(i64==4)
                        {
                              if(QLinks<0.5)
                                    QLinks      = QLinks * (1.67 - QLinks * 1.078);
                              else
                                    QLinks      = QLinks * 0.869 + 0.131;
                        }
                  }
                  break;
            }
      }


    //---QRechts-ermitteln---------------------------------
      x     = 9;
      y     = 1;

      for(i=0; i<100; i++)
      {
            ix    = BRM_nint(x) - 1;
            iy    = BRM_nint(y) - 1;

            for(j=0; j<i64; j++)
            {
                  a     = BRM_g[g64][j][ix][iy];
                  x     += s[j] * a;
                  y     += c[j] * a;
            }

            if(x>9)
            {
                  ARecht      = 0;
                  QRecht      = 0;
                  break;
            }

            if( x<1.01 || y<1 )
            {
                  ARecht      = 1;
                  QRecht      = 1;
                  break;
            }

            if(y>8.95)
            {
                  if(x>8.98)
                  {
                        ARecht      = 0;
                        QRecht      = 0;
                  }
                  else
                  {
                        ARecht      = 1;
                        QRecht      = 1 - (x - 1) / 8;

                        if(i64==4)                          
                        {
                              if(QRecht<0.5)
                                    QRecht      = QRecht * (1.67 - QRecht * 1.078);
                              else
                                    QRecht      = QRecht * 0.869 + 0.131;
                        }
                  }
                  break;
            }
      }

      //---QMitte-ist-Rest-/-QBinaer'Bits-setzen-------------
      QMitte      = 1 - QLinks - QRecht;
      QBinaer     = ALinks + AMitte + ARecht;

      if(QMitte<=0.01)
      {
            a           = QLinks + QRecht;
            QLinks      /= a;
            QRecht      /= a;
            QMitte      = 0;
            QBinaer     -= 2;
      }
}

Generated by  Doxygen 1.6.0   Back to index