revert using :third_order for calc_z

Author: orso82

src/math.jl

@@ -740,8 +740,9 @@ The finite difference `method` of the gradient can be one of [:third_order, :sec
 
 NOTE: the inverse scale length is NEGATIVE for typical density/temperature profiles
 """
-function calc_z(x::AbstractVector{<:Real}, f::AbstractVector{<:Real}; method::Symbol=:backward)
-    return gradient(x, f; method) ./ f
+function calc_z(x::AbstractVector{<:Real}, f::AbstractVector{<:Real}; method::Symbol=:third_order)
+    g = gradient(x, f; method)
+    return g ./ f
 end
 
 """
@@ -750,16 +751,15 @@ end
 Backward integration of inverse scale length vector with given edge boundary condition
 """
 function integ_z(rho::AbstractVector{<:Real}, z_profile::AbstractVector{<:Real}, bc::Real)
-    f = interp1d(rho, z_profile, :cubic)
     profile_new = similar(rho)
     profile_new[end] = bc
     for i in length(rho)-1:-1:1
         a = rho[i]
         fa = z_profile[i]
         b = rho[i+1]
         fb = z_profile[i+1]
-        simpson_integral = (b - a) / 6 * (fa + 4 * f((a + b) * 0.5) + fb)
-        profile_new[i] = profile_new[i+1] * exp(simpson_integral)
+        trapz_integral = (b - a) * (fa + fb) / 2.0
+        profile_new[i] = profile_new[i+1] * exp(trapz_integral)
     end
     return profile_new
 end
@@ -1060,7 +1060,7 @@ function perimeter(r::AbstractVector{T}, z::AbstractVector{T})::T where {T<:Real
     end
 
     # If open, add distance from last point to first point
-    if is_open_polygon(r,z)
+    if is_open_polygon(r, z)
         dx = r[1] - r[end]
         dy = z[1] - z[end]
         perimeter += sqrt(dx^2 + dy^2)

src/physics/transport.jl

@@ -8,7 +8,8 @@
 
 Updates profile_old with the scale lengths given by z_transport_grid
 
-If `rho_ped` > transport_grid[end]` then scale-length is linearly interpolated between `transport_grid[end]` and `rho_ped`
+If `rho_ped > transport_grid[end]` then scale-length is linearly interpolated between `transport_grid[end]` and `rho_ped`
+if `rho_ped < transport_grid[end]` then scale-length then boundary condition is at `transport_grid[end]`
 """
 function profile_from_z_transport(
     profile_old::AbstractVector{<:Real},
@@ -17,24 +18,24 @@ function profile_from_z_transport(
     z_transport_grid::AbstractVector{<:Real},
     rho_ped::Real=0.0)
 
-    z = -calc_z(rho, profile_old)
+    z = calc_z(rho, profile_old)
 
-    transport_idices = [argmin(abs.(rho .- rho_x)) for rho_x in transport_grid]
-    index_nml = argmin(abs.(rho .- rho_ped))
-    if index_nml > transport_idices[end]
-        transport_idices = vcat(1, transport_idices, index_nml)
-        z_transport_grid = vcat(0.0, z_transport_grid, -z[index_nml])
+    transport_indices = [argmin(abs.(rho .- rho_x)) for rho_x in transport_grid]
+    index_ped = argmin(abs.(rho .- rho_ped))
+    index_last = transport_indices[end]
+    if index_ped > index_last
+        transport_indices = vcat(1, transport_indices, index_ped)
+        z_transport_grid = vcat(0.0, z_transport_grid, z[index_ped])
     else
-        transport_idices = vcat(1, transport_idices)
+        transport_indices = vcat(1, transport_indices)
         z_transport_grid = vcat(0.0, z_transport_grid)
     end
-    rho_transport_grid = rho[transport_idices]
 
-    z[1:transport_idices[end]] = -interp1d(rho_transport_grid, z_transport_grid).(rho[1:transport_idices[end]])
+    z[1:index_last] = interp1d(transport_indices, z_transport_grid).(1:index_last)
 
     profile_new = similar(profile_old)
-    profile_new[transport_idices[end]:end] = profile_old[transport_idices[end]:end]
-    profile_new[1:transport_idices[end]] = integ_z(rho[1:transport_idices[end]], z[1:transport_idices[end]], profile_new[transport_idices[end]])
+    profile_new[index_last:end] = @views profile_old[index_last:end]
+    profile_new[1:index_last] = @views integ_z(rho[1:index_last], -z[1:index_last], profile_new[index_last])
 
     return profile_new
 end