anchor.diff:          turbI = sorted_x_idx(d) + 1
           
           do, p=1, nRotorPoints
  +        
  +            ! initialize the TI_area_ratio to 0.0 for each turbine
  +            ti_area_ratio = 0.0_dp
       
               ! initialize deficit summation term to zero
               deficit_sum = 0.0_dp
  
positive.code:     
          ! final velocity calculation for turbine turbI (average equally across all points)
          wtVelocity(turbI) = wtVelocity(turbI)/nRotorPoints
         
      end do
  
      !print *, "here 0"
      
      !print *, "entering visualization point loops"
      
      do, loc=1, nSamples
           
          ! set combined velocity deficit to zero for loc
          deficit_sum = 0.0_dp
          
          !print *, "entering visualization point turbine loop"
          
          do, turb=1, nTurbines ! at turbineX-locations
              !print *, "loc = ", loc, "turb = ", turb
          
              ! set velocity deficit at loc due to turb to zero
              deltav = 0.0_dp
  
              ! determine the onset location of far wake
              call x0_func(rotorDiameter(turb), yaw(turb), Ct(turb), alpha, TIturbs(turb), &
                          & beta, x0)
              !print *, "x0: ", x0/rotorDiameter(turb)
          
              ! determine the initial wake angle at the onset of far wake
              call theta_c_0_func(yaw(turb), Ct(turb), theta_c_0)
          
              ! downstream distance between turbines
              x = velX(loc) - turbineXw(turb)
                  
              ! downstream distance from far wake onset to downstream turbine
              deltax0 = x - x0
              
              !print *, "entering logic statements"
              
              if (x > 0.0_dp) then
              
                  !print *, "entering first if logic"
              
                  ! horizontal spread
                  call sigmay_func(ky, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
              
                  ! vertical spread
                  call sigmaz_func(kz, deltax0, rotorDiameter(turb), sigmaz)
              
                  ! horizontal cross-wind wake displacement from hub
                  call wake_offset_func(rotorDiameter(turb), theta_c_0, x0, yaw(turb), &
  
negative.code:             
                  ! vertical distance from downstream hub location to wake center
                  deltaz = velZ(loc) - turbineZ(turb)
                                  
                  ! find the final point where the original model is undefined
                  call discontinuity_point_func(x0, rotorDiameter(turb), ky_local(turb), kz_local(turb), yaw(turb), Ct(turb), & 
                                               & discontinuity_point)
  