diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 979cf0d490fd07f811631b1829892e868e67bb00..3898c227cbbad19d6603d6fcea992de357ed02fc 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -190,7 +190,7 @@ function L1L2TVState{m}(mesh; T, tdata, S,
 end
 
 function OptFlowState(mesh;
-	alpha1, alpha2, beta, lambda, gamma1, gamma2)
+	alpha1, alpha2, beta, lambda, gamma1, gamma2, Schoice)
     alpha2 > 0 || beta > 0 ||
         throw(ArgumentError("operator B is singular with these parameters"))
 
@@ -208,7 +208,7 @@ function OptFlowState(mesh;
     # tdata will be something like nabla(fw)
     T(tdata, u) = tdata * u
     T(::typeof(adjoint), tdata, v) = tdata' * v
-    S(u, nablau) = nablau'
+    S(u, nablau) = Schoice == :nabla ? nablau' : u
 
     est = FeFunction(Vest, name = "est")
     g = FeFunction(Vg, name = "g")
@@ -1526,7 +1526,7 @@ function optflow(ctx)
     st = OptFlowState(mesh;
         ctx.params.alpha1, ctx.params.alpha2,
         ctx.params.lambda, ctx.params.beta,
-        ctx.params.gamma1, ctx.params.gamma2)
+        ctx.params.gamma1, ctx.params.gamma2, ctx.params.Schoice)
 
     function warp!()
         println("warp and reproject ...")
@@ -1582,8 +1582,8 @@ function optflow(ctx)
             println("norm_step = $norm_step_")
 
             # interior newton stop criterion
-            norm_step_ > eps_newton && k_newton < 30 && continue
-            k_newton >= 30 && @warn "Newton reached maximum number of iterations"
+            norm_step_ > eps_newton && k_newton < ctx.params.newton_max_iters && continue
+            k_newton >= ctx.params.newton_max_iters && @warn "Newton reached maximum number of iterations"
             ctx.params.warp || break
             k_newton = 0
 
@@ -1642,7 +1642,7 @@ function optflow(ctx)
     saveimg(joinpath(ctx.outdir, "fw.png"), to_img(imgfw))
     savedata(joinpath(ctx.outdir, "data.tex");
         eps_newton, eps_warp, ctx.params.n_refine,
-        st.alpha1, st.alpha2, st.lambda, st.beta, st.gamma1, st.gamma2,
+        st.alpha1, st.alpha2, st.lambda, st.beta, st.gamma1, st.gamma2, ctx.params.Schoice,
         width=size(u_sampled, 2), height=size(u_sampled, 3),
         endpoint_error_mean, endpoint_error_stddev,
         angular_error_mean, angular_error_stddev,
@@ -1667,8 +1667,7 @@ function experiment_optflow_middlebury_all_benchmarks(ctx)
         #example == "Dimetrodon" && continue
         ctx(experiment_optflow_middlebury, example;
             alpha1 = 10., alpha2 = 0., lambda = 1., beta = 1e-5,
-            gamma1 = 1e-4, gamma2 = 1e-4)
-        #return
+            gamma1 = 1e-4, gamma2 = 1e-4, Schoice = :nabla, newton_max_iters = 30)
     end
 end
 
@@ -1694,6 +1693,41 @@ function experiment_optflow_middlebury_warping_comparison_adaptive(ctx)
         warp = true, refine = true, n_refine = 6)
 end
 
+function experiment_optflow_schoice(ctx)
+    function run(beta, Schoice)
+        example = "Dimetrodon"
+        duration = @elapsed begin
+            st = ctx(experiment_optflow_middlebury, example;
+                alpha1 = 10., alpha2 = 0., lambda = 1., beta,
+                gamma1 = 1e-4, gamma2 = 1e-4, Schoice,
+                warp = true, refine = false, n_refine = 0, newton_max_iters = 15)
+        end
+
+        u_sampled = sample(st.u)
+        u_flow = to_img(u_sampled)
+
+        # not yet in ctx
+        gtflow = FileIO.load(joinpath(ctx.indir, example, "flow10.flo"))
+        endpoint_errors = collect(skipmissing(map_flow(endpoint_error, u_flow, gtflow)))
+        angular_errors = collect(skipmissing(map_flow(angular_error, u_flow, gtflow)))
+        endpoint_error_mean = mean(endpoint_errors)
+        endpoint_error_stddev = stdm(endpoint_errors, endpoint_error_mean)
+        angular_error_mean = mean(angular_errors)
+        angular_error_stddev = stdm(angular_errors, angular_error_mean)
+
+        return (;Schoice, beta, duration, endpoint_error_mean, endpoint_error_stddev, angular_error_mean, angular_error_stddev)
+    end
+
+    df = DataFrame()
+    for Schoice in (:id, :nabla)
+        for beta in 10 .^ range(-6,3,10)
+            res = run(beta, Schoice)
+            push!(df, res)
+        end
+    end
+    CSV.write(joinpath(ctx.outdir, "schoice.csv"), df)
+end
+
 
 function test_image(n = 2^6; supersample_factor = 16)
     q = supersample_factor