anchor.diff:  end subroutine porteagel_analyze
   
   ! implementation of the Bastankhah and Porte Agel (BPA) wake model for visualization
  -subroutine porteagel_visualize(nTurbines, nSamples, nRotorPoints, turbineXw, sorted_x_idx, turbineYw, turbineZ, &
  +subroutine porteagel_visualize(nTurbines, nSamples, nRotorPoints, nCtPoints, turbineXw, &
  +                             sorted_x_idx, turbineYw, turbineZ, &
                                rotorDiameter, Ct, wind_speed, &
                                yawDeg, ky, kz, alpha, beta, TI, RotorPointsY, RotorPointsZ, & 
                                z_ref, z_0, shear_exp, velX, velY, velZ, &
                                wake_combination_method, TI_calculation_method, &
  -                             calc_k_star, opt_exp_fac, wsArray)
  +                             calc_k_star, opt_exp_fac, wake_model_version, &
  +                             use_ct_curve, ct_curve_wind_speed, ct_curve_ct, wsArray)
                                
       ! independent variables: turbineXw turbineYw turbineZ rotorDiameter
       !                        Ct yawDeg
  
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:                     if (x > discontinuity_point) then
      
                          ! velocity difference in the wake
                          call deltav_func(deltay, deltaz, Ct(turb), yaw(turb), &
                                          & sigmay, sigmaz, rotorDiameter(turb), & 
                                          & wake_model_version, kz_local, x, &
                                          & opt_exp_fac, deltav)  
                                          
                      ! near wake region (linearized)
                      else
                          
  