diff --git a/docs/src/tutorials.md b/docs/src/tutorials.md
index cba8128..3eab39f 100644
--- a/docs/src/tutorials.md
+++ b/docs/src/tutorials.md
@@ -271,6 +271,340 @@ end
 
 Given the finite slippage we do not observe *large* deviations from the $\alpha=1/10$ powerlaw in the long time limit.
 
+
+## Active film
+
+In [Richter et all](https://arxiv.org/abs/2402.14635) we introduced a thin film model for chemical reactions inside a liquid film aimed at Supported liquid phase catalysts. For the details of the model look at the article. Here we give a minimal example of how to run a simulation, starting from a homogenous film, leading to film rupture. 
+
+Lets get straight in it 
+```julia 
+using Swalbe, DelimitedFiles
+function simulation(;
+    D=5e-6,			# Diffusion coefficient
+    Tmax=1e3, 			# Simulation time
+    tdump=Int(floor(Tmax/100)), # dump time
+    h= 1,  			# initial film height
+    rho=1,			# initial catalyst concentration
+    ϵ=0.001,			# initial noise
+    γ_0=0.001,			# reference surface tension in the absence of chemical compounds
+    Γ = 0.1, 			# Surface tension effect product rho_A
+    GammaB=0.1,			# Surface tension effect reactant rhoB
+    data = "data",		# location to save data
+    L=2^12, 			# system size
+    alpha=0.0, 			# catalyst vertical distribution parameter 
+    theta_0=1/9,		# reference contact angle in fractions of pi
+    n=9, m=3, hmin=0.1, hcrit=0.05, hcrit2=0.05,	# disjoining pressure parameters
+    delta=1.0,			# slip length
+    production_rate=0.1,	# production rate
+    sigma_A_up=0.1, 		# evaporation rate product rho_A
+    sigma_B_up=0.1, 		# evaporation rate reactant rho_B
+    rho_BRes_sigma_B_down=0.1,	# sorption rate reactant rho_B
+    rho_0B=rho_BRes_sigma_B_down*h/(production_rate*rho + sigma_B_up),	# initial reactant concentration 
+    rho_0A=production_rate*rho*rho_0B/sigma_A_up,			# initial product cocentration 
+    D_Product=D*20, D_B=D*20,	# Diffusion coefficients Product rho_A and reactant rho_B
+    rho_crit=0.05,		# precursor length catalyst
+    mu=1/6,			# viscosity
+)
+	# Set up sys consts
+    	    GammafacA=rho_0A/(h*γ_0)
+    	    GammafacB=rho_0B/(h*γ_0)
+	    sys = Swalbe.SysConstActive_1D{Float64}(
+        	L =  L,
+        	D=D,
+        	Γ = Γ/GammafacA,
+        	GammaB=GammaB/GammafacB,
+        	Tmax=Tmax,
+        	tdump=tdump,
+        	γ_0 = γ_0,
+        	δ=delta,
+        	alpha=alpha,
+        	θ_0 =theta_0,
+        	n=n,m=m,hmin=hmin,hcrit=hcrit, hcrit2=hcrit2,
+        	production_rate=production_rate,
+        	sigma_A_up=sigma_A_up,
+        	sigma_B_up=sigma_B_up,
+        	rho_BRes_sigma_B_down=rho_BRes_sigma_B_down,
+        	D_Product=D_Product,
+        	D_B=D_B,
+        	rho_crit=rho_crit,
+        	b_A=rho_crit,
+        	b_B=rho_crit,
+        	μ=mu
+	)
+	# set up simulation fields 
+	state = Swalbe.Sys(sys)
+	# Set initial values
+	state.rho.=rho
+        state.rho_A .= rho_0A
+        state.rho_B .= rho_0B
+	# This set up the film height as constant plus a small noise without high wavenumber frequencies
+        Swalbe.browniannoise!(state.height, h, ϵ, L/2-L/10)
+	i = 0
+	# The LBM time loop
+	for t in 0:sys.Tmax
+		# dump data
+		if t % sys.tdump == 0 
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_height.csv", state.height)
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_rho.csv", state.rho)
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_rho_A.csv", state.rho_A)
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_rho_B.csv", state.rho_B)
+			i += 1
+			println("Time step t=$t")
+		end
+		# system update, mostly as in a standard TFE simulation, the parts that are different are commented
+		Swalbe.surface_tension_4_fields!(state, sys)
+		Swalbe.filmpressure_fast!(state, sys)
+		Swalbe.h∇p!(state)
+            	Swalbe.slippage!(state, sys)
+		# claculate the Marangoni stress
+		Swalbe.∇γ!(state, sys)
+		# Include the Marangoni stress into the forcing
+            	state.F .= -state.h∇p .- state.slip .+ state.stress
+		# Do the LBM loops for catalys state.rho, product state.rho_A, and reactant state.rho_B
+		Swalbe.update_rho_LBM!(state,sys)
+		Swalbe.update_rho_A_LBM_precursor!(state,sys)
+		Swalbe.update_rho_B_LBM_precursor!(state,sys)
+		Swalbe.equilibrium!(state)
+		Swalbe.BGKandStream!(state, sys)
+		Swalbe.moments!(state)	
+	end
+end
+
+# run the simulation 
+simulation()
+```
+
+
+That's it
+
+## Thin films on curved substrate
+
+It is possible to include a curved substrate into the model using the substrate profile `b` into the pressure
+```math
+p_{\text{cap}} = -(1 + 0.5 \nabla s^2) \gamma (1 - \cos(\theta)) \phi'(h) - \gamma \nabla^2 (h+b)
+```
+
+Here is an example on how to do this 
+
+```julia 
+using Swalbe, DelimitedFiles, Dates
+function simulation(;
+	Tmax=1e8, 			# Simulation time
+	tdump=Int(floor(Tmax/100)),	# dump time
+	h= 1,				# initial film height  
+	eps=0.001, 			# initial noise
+	gamma=1/6^2,			# surface tension
+	data = "data",			# folder for saving data
+ 	L=2^12,				# system size
+    	theta=1/9,			# contact angle
+    	n=9, m=3, hmin=0.1, 		# parameters for disjoining pressure
+	delta=1.0,			# slip length, in combination with heigher hmin i
+	omega=3,			# number of periods of the underlying substrate
+)
+	# set up system
+    	sys = Swalbe.SysConst_1D{Float64}(L =  L ,Tmax=Tmax, tdump=tdump, γ = gamma, δ=delta,θ =theta, n=n,m=m,hmin=hmin)
+     	try
+        	mkdir(data)
+            	mkdir("$(data)/figs")
+        catch
+        end
+    	state = Swalbe.Sys(sys, kind="curved")
+        # Initial data
+    	state.height.= h .+ eps .* rand(sys.L)
+	# set up system profile
+    	state.substrate .=[1/4*sin(i*omega*2*pi/sys.L) for i in 1:sys.L]
+	# We need the gradient of the substrate profile
+    	Swalbe.∇f!(state.grad_substrate, state.substrate, state.dgrad)
+    	println("Starting Lattice Boltzmann time loop for active thin film")
+    	i=0
+    	for t in 0:sys.Tmax
+		# save data
+        	if t % sys.tdump == 0
+                	writedlm("$(data)/$(Swalbe.to_4_digits(i))_height.csv", state.height)
+            		i+=1
+                end
+		# system update, using the pressure for curved substrates
+        	Swalbe.filmpressure_curved!(state, sys)
+        	Swalbe.h∇p!(state)
+        	Swalbe.slippage!(state, sys)
+        	state.F .= -state.h∇p .- state.slip
+        	Swalbe.equilibrium!(state)
+        	Swalbe.BGKandStream!(state, sys)
+        	Swalbe.moments!(state)
+    	end
+    println("Finished timeloop")
+    return nothing
+end
+simulation()
+```
+
+
+## Miscible liquids
+
+Here I show to simulate the colaescence of two droplets of miscible liquids with different surface tensions, lets say water and ethanol. 
+
+```julia
+using Swalbe, DelimitedFiles,  Dates
+
+
+function run_Miscible(;
+    eps=0.01,                                   # initial noise
+    Tmax=Int(1e7),                              # Simulation time
+    dumps=100,                                  # number of dumps
+    tdump=Int(floor(Tmax/dumps)),               # dump time
+    L=2^9,                                      # system size
+    delta=1.0,                                  # slip length
+    hmin=0.1, hcrit=0.05, n=9, m=3,             # disjoining pressure parameters
+    gamma= 0.01 .* ones(3,2),                   # surface tensions
+    data="path/to/data",                        # save path
+    D=1e-2,                                     # diffusion coeficient
+    r=50,                                       # Initial droplet radiusj
+    center_weight=0.5                           # weight for surface tension smoothing
+    )
+    #Make folders
+    datathis= "$(data)/Tmax=$Tmax,L=$L,delta=$delta,gamma=$(replace(replace("$gamma", ";"=>","), " " => "_")),D=$D,r=$r,hmin=$hmin,hcrit=$hcrit"
+        try
+        mkdir(dataall)
+        catch
+        end
+        try
+        mkdir(datathis)
+        catch
+        end
+        try
+            mkdir("$(datathis)/figs")
+        catch
+        end
+    #Setup system
+    sys = Swalbe.SysConstMiscible_1D(L=L, Tmax=Tmax, tdump=tdump, gamma=gamma, delta=delta, n=n, m=m, D=D, hmin=hmin, hcrit=hcrit)
+    state = Swalbe.Sys(sys)
+    # Initial data
+    theta_1=acos((sys.gamma[3,1]-sys.gamma[2,1])/sys.gamma[1,1])
+    theta_2=acos((sys.gamma[3,2]-sys.gamma[2,2])/sys.gamma[1,2])
+    state.height[:,1] .= Swalbe.droplet_base(state.height[:,1], r, theta_1/pi, Int(sys.L/2 - r/2), precursor=(sys.hmin-2*sys.hcrit)/2)
+    state.height[:,2] .= Swalbe.droplet_base(state.height[:,2], r, theta_2/pi, Int(sys.L/2 + r/2), precursor=(sys.hmin-2*sys.hcrit)/2)
+    bridge_height=[]
+    i=0
+    println("Starting LBM time loop")
+    for t in 0:sys.Tmax
+        if t % sys.tdump == 0
+            #Store data
+            writedlm("$(datathis)/$(Swalbe.to_4_digits(i))_h_t=$t",state.height)
+            println("time step t=$t")
+            i+=1
+        end
+        #actual simulations
+            # surface tension calculation from the height fields
+            # Swalbe.surface_tension!(state,sys)
+            # surface tension calculation from the height fields with moving averadge smoothing applied
+            Swalbe.surface_tension_smooth!(state,sys, center_weight=center_weight)
+            Swalbe.filmpressure!(state, sys)
+            Swalbe.h∇p!(state,sys)
+            Swalbe.slippage!(state,sys)
+            # calculate Marangoni shear stress
+            Swalbe.∇gamma!(state, sys)
+            state.F .= -state.h∇p  .-  state.slip .+ state.stress
+            Swalbe.equilibrium!(state)
+            Swalbe.BGKandStream!(state, sys)
+            Swalbe.moments!(state, sys)
+    end
+end
+
+
+dataall = "data/$(Dates.today())Miscible_coal_test_00"
+gamma=ones(3,2)
+gamma_0=0.002
+# fixed contact angles
+Delta=0.2
+gamma[1,1]=gamma_0*(1+Delta/2)
+gamma[1,2]=gamma_0*(1-Delta/2)
+# solid vapour
+gamma[3,1]=2*gamma_0
+gamma[3,2]=2*gamma_0
+theta=1/9
+# fixed contact angle
+gamma[2,1]=gamma[3,1]-gamma[1,1]*cospi(theta)
+gamma[2,2]=gamma[3,2]-gamma[1,2]*cospi(theta)
+run_Miscible(Tmax=Int(1e7),gamma=gamma, D=1.5e-4, r=200, delta=1.0, hmin=0.2, L=2^10, data=dataall)
+```
+
+Have fun with it
+
+## Multilayer thin films
+
+This shows how to simulate a layer film resting on a layer of water. This we call a two layer multilayer system. To read more about the theory see [Richter et al.](https://arxiv.org/abs/2409.16659). 
+
+This is the code on how to do this: 
+```julia 
+using Pkg
+
+Pkg.activate(".")
+
+
+using Swalbe, DelimitedFiles, CUDA
+
+
+function run_multilayer(;
+    h_1=2,		# initial height lower layer
+    h_2=1,		# initial height upper layer
+    eps=[0.001,0.001],	# initial noise
+    Tmax=Int(1e7),	# simulation time
+    dumps=100,		# number of dumps
+    tdump=Int(floor(Tmax/dumps)), 	# dump time
+    Lx=2^8, Ly=2^8,	# system size
+    delta=1.0,		# slip lower layer
+    delta_2=1.0,	# slip upper layer
+    n=9, m=3,		# disjoining pressure exponents
+    gamma= [0 0.001 0.001 0.001 ; 0 0 0.001 0.001 ; 0 0 0 0.001; 0 0 0 0],	# surface tensions gamma[i,j] is the interaction between layer i-1 and j-1 having 0 for the solid substrat and 3 the atmosphere
+    data="data", 	# path to save simulation dumps
+    hmin=0.1,		# precursor layer
+    mu=[1/6 1/6],	# viscosities 
+    )
+    #Setup system
+    sys = Swalbe.SysConstMultiLayer(Lx=Lx, Ly=Ly, Tmax=Tmax, tdump=tdump, gamma=gamma, delta=[delta,delta_2], n=n, m=m, hmin=hmin, mu=mu)
+    state = Swalbe.Sys(sys, device=device)
+    # this is needed to send back and force in between GPU and CPU
+    host_height=zeros(Lx,Ly,2)
+    #initial conditions
+    host_height[:,:,1] .= Swalbe.browniannoise2D!(host_height[:,:,1], h_1, eps[1], Lx/2-Lx/10)
+    host_height[:,:,2] .= Swalbe.browniannoise2D!(host_height[:,:,2], h_2, eps[2], Lx/2-Lx/10)
+    CUDA.copyto!(state.height, host_height)
+    i=0
+    println("Starting LBM time loop")
+    for t in 0:sys.Tmax
+        if t % sys.tdump == 0 
+            #Store data
+                CUDA.copyto!(host_height,state.height)
+                writedlm("$(datathis)/$(Swalbe.to_4_digits(i))_h1_t=$t",host_height[:,:,1])
+                writedlm("$(datathis)/$(Swalbe.to_4_digits(i))_h2_t=$t",host_height[:,:,2])
+            #status report
+            println("Time step t=$t)
+            i+=1
+        end
+        #actual simulations
+        if run 
+            Swalbe.filmpressure!(state, sys)
+	    Swalbe.h∇p!(state,sys)
+            Swalbe.slippage!(state, sys)
+	    state.Fx .= -state.h∇px  .+  state.slipx
+	    state.Fy .= -state.h∇py  .+  state.slipy
+            Swalbe.equilibrium!(state)
+            Swalbe.BGKandStream!(state)
+            Swalbe.moments!(state, sys)
+        end
+    end 
+end
+
+
+dataall = "path\to\save\dumps"
+gamma_0=0.01
+gamma=gamma_0*ones(4,4)
+gamma[2,4]=gamma_0*(1+cospi(2/9))
+gamma[1,3]=gamma_0*(2.1)
+gamma[1,4]=gamma[2,4]+gamma[1,3]-gamma[2,3]
+run_multilayer(Tmax=Int(1e6),gamma=gamma, h_2=1.0, h_1=6.0, delta=1.0, delta_2=1.0, Lx=2^8,Ly=2^8, mu=[0.1 0.1], hmin=0.2, devicenumber=2)
+```
+
 ## Further tutorials
 
 More tutorials will follow in the future.
@@ -279,4 +613,4 @@ So be sure to check out the docs every now and then.
 
 The next tutorial will be about switchable substrates. 
 In this case the wettability can not only addressed locally but also with a time dependency.
-Here is what happens if the time frequency is [high](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) and this happens if we update with a [lower](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) frequency.
\ No newline at end of file
+Here is what happens if the time frequency is [high](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) and this happens if we update with a [lower](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) frequency.
diff --git a/scripts/init.jl b/scripts/init.jl
new file mode 100644
index 0000000..e664e9e
--- /dev/null
+++ b/scripts/init.jl
@@ -0,0 +1,5 @@
+using Pkg
+
+Pkg.activate(".")
+
+using Swalbe
diff --git a/src/differences.jl b/src/differences.jl
index 5425344..5993c4f 100644
--- a/src/differences.jl
+++ b/src/differences.jl
@@ -229,6 +229,20 @@ function ∇f!(output, f::Vector, dgrad)
     return nothing
 end
 
+
+function ∇f_Miscible!(output, f::Array, dgrad)
+    fip, fim = viewneighborsMiscible_1D(dgrad)
+    # One dim case, central differences
+    circshift!(fip, f, 1)
+    circshift!(fim, f, -1)
+    # In the end it is just a weighted sum...
+    output .= -0.5 .* (fip .- fim)
+    return nothing
+end
+
+
+
+
 """
     viewneighbors(f)
 
diff --git a/src/forcing.jl b/src/forcing.jl
index e9cc663..08ddca5 100644
--- a/src/forcing.jl
+++ b/src/forcing.jl
@@ -80,8 +80,8 @@ slippage!(state::Expanded_1D, sys::Consts_1D) = slippage!(
     sys.param.μ,
 )
 
-slippage!(state::Active_1D, sys::SysConstsActive_1D) = slippage!(state.slip, state.height, state.vel, sys.δ, sys.μ)
-slippage!(state::LBM_state_1D, sys::SysConstsMiscible_1D) = slippage!(state.slip, state.height, state.vel, sys.delta, sys.μ)
+slippage!(state::Active_1D, sys::SysConstActive_1D) = slippage!(state.slip, state.height, state.vel, sys.δ, sys.μ)
+slippage!(state::LBM_state_1D, sys::SysConstMiscible_1D) = slippage!(state.slip, state.height, state.vel, sys.delta, sys.μ)
 
 # Dirty hack for reducing slip length
 function slippage2!(state::LBM_state_2D, sys::SysConst)
@@ -1579,3 +1579,30 @@ function view_four(f)
 
     return f1, f2, f3, f4
 end
+
+
+"""
+    surface_tension_gradient!(state)
+
+Computes the gradient of a spatially resolved surface tension field via central differneces and then concludes the surface stress as
+
+`` \\frac{h^2 + 2 \\delta h }{2M(h)}``
+
+with
+
+``M(h)= \\frac{2h^2 + 6 \\delta h + 3\\delta^2 }{ 6}``
+
+and implicitly ``\\rho_h=1``.
+"""
+function ∇γ!(state::Active_1D, sys::SysConstActive_1D)
+    fip, fim = viewneighbors_1D(state.dgrad)
+    # One dim case, central differences
+    circshift!(fip, state.γ, 1)
+    circshift!(fim, state.γ, -1)
+    # In the end it is just a weighted sum...
+    # state.∇γ .= -3/2 .* ((fip .- fim) ./ 2.0)
+    state.∇γ .=((fip .- fim) .* -0.5)
+    state.stress .= state.∇γ .* (state.height .* state.height .* 0.5 .+ sys.δ .* state.height) .* 6 ./ (2 .* state.height .* state.height .+ 6*sys.δ .* state.height .+ 3*sys.δ*sys.δ )
+    return nothing
+end
+
diff --git a/src/initialize.jl b/src/initialize.jl
index 46fafe1..eee43ea 100644
--- a/src/initialize.jl
+++ b/src/initialize.jl
@@ -4,11 +4,13 @@ abstract type LBM_state end
 abstract type LBM_state_2D <: LBM_state end
 
 abstract type Expanded_2D <: LBM_state_2D end
+abstract type MultiLayer_2D <: LBM_state_2D end
 # Derived type for 1D LBM states
 abstract type LBM_state_1D <: LBM_state end
 
 abstract type Expanded_1D <: LBM_state_1D end
 abstract type Boundary_1D <: Expanded_1D end
+abstract type Active_1D <: LBM_state_1D end
 
 # Parent type for system specific constants
 abstract type Consts end
@@ -1238,7 +1240,6 @@ function Sys(sysc::Consts_1D; T = Float64, kind = "simple")
             ∇γ = zeros(sysc.L),
             fbound = zeros(sysc.L, 3),
         )
-    end
     elseif kind == "curved"
             if device=="CPU"
             	    dyn = State_curved_1D{T}(
diff --git a/src/initialvalues.jl b/src/initialvalues.jl
index 18fe942..d9f5884 100644
--- a/src/initialvalues.jl
+++ b/src/initialvalues.jl
@@ -445,3 +445,14 @@ function trianglepattern(
     return θ, P
 end
 
+function droplet_base(height::Vector, s, theta, center; precursor=0.05)
+	L=size(height)[1]
+	dummy=zeros(L)
+	r=s/(2*sinpi(theta))
+	h=r*(1-cospi(theta))
+	height .= precursor
+	for x in Int(center-s/2):Int(center+s/2)
+		dummy[x]=sqrt(r^2-(x-center)^2)-r+h
+	end
+	return dummy.=max.(dummy, precursor)
+end
diff --git a/src/logic.jl b/src/logic.jl
index 8ff67db..8ee7a9f 100644
--- a/src/logic.jl
+++ b/src/logic.jl
@@ -260,7 +260,7 @@ end
 
 
 """
-    fast_disj_32(arg::Float64)
+    fast_32(arg::Float64)
 
 Computes `arg^3 - arg^2` more efficiently using the `power_3` and `power_2` functions. This method avoids computing `arg^2` multiple times.
 
@@ -272,22 +272,18 @@ Computes `arg^3 - arg^2` more efficiently using the `power_3` and `power_2` func
 
 # Example
 ```jldoctest
-julia> fast_disj_32(2)
+julia> fast_32(2)
 4
 ```
 """
-function fast_disj_32(arg::Float64)
+function fast_32(arg::Float64)
     return power_3(arg) - power_2(arg)
 end
 
-function fast_disj_32(arg::Float64, fac::Float64)
-    # This was introducing anisotropy for reasons I don't understand
-    # return fac*power_3(arg)-power_2(arg)
-    return power_3(arg)-power_2(arg)
-end
+
 
 """
-    fast_disj_93(arg::Float64)
+    fast_93(arg::Float64)
 
 Computes `arg^3 - arg^3(arg)` more efficiently using the `power_3` function twice. This method avoids redundant power calculations.
 
@@ -299,21 +295,15 @@ Computes `arg^3 - arg^3(arg)` more efficiently using the `power_3` function twic
 
 # Example
 ```jldoctest
-julia> fast_disj_93(2)
+julia> fast_93(2)
 56
 ```
 """
-function fast_disj_93(arg::Float64)
+function fast_93(arg::Float64)
     temp = power_3(arg)
     return power_3(temp) - temp
 end
 
-function fast_disj_93(arg::Float64, fac::Float64)
-    temp = power_3(arg)
-    # # This was introding anisotropy for reasons I don't understand
-    # return fac*power_3(temp)-(fac+1)/2*temp
-    return power_3(temp)-temp
-end
 
 """
     to_4_digits(i)
@@ -350,7 +340,7 @@ function to_4_digits(i)
     end
 end
 
-function fast_disj_63(arg::Float64, fac::Float64)
+function fast_63(arg::Float64, fac::Float64)
     temp = power_3(arg)
     # This was introducing anisotropy for reasons I don't understand
     # return fac*power_2(temp)-temp
diff --git a/src/pressure.jl b/src/pressure.jl
index 8b5ebab..0aeabb7 100644
--- a/src/pressure.jl
+++ b/src/pressure.jl
@@ -392,15 +392,15 @@ function filmpressure!(state::MultiLayer_2D, sys::SysConstMultiLayer)
 	if sys.n==9 && sys.m==3
        # if sys.layers == 2
        # store the one field we need twice
-       state.hi[:,:,1].=(sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_disj_93.(sys.hmin./(state.height[:,:,2] .+ sys.hcrit)).* (1 .+ state.grad_h_sq ./2)
+       state.hi[:,:,1].=(sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_93.(sys.hmin./(state.height[:,:,2] .+ sys.hcrit)).* (1 .+ state.grad_h_sq ./2)
         # disjoining pressure
         state.pressure[:,:,1] .=  (
-                (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])                  *(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_disj_93.(sys.hmin./(state.height[:,:,1] .+ sys.hcrit))
+                (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])                  *(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_93.(sys.hmin./(state.height[:,:,1] .+ sys.hcrit))
                 .- state.hi[:,:,1]
             )
             state.pressure[:,:,2] .=  (
                 state.hi[:,:,1]
-                .+ (sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_disj_93.(2*sys.hmin./(state.height[:,:,1] .+ state.height[:,:,2] .+ 2*sys.hcrit))
+                .+ (sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_93.(2*sys.hmin./(state.height[:,:,1] .+ state.height[:,:,2] .+ 2*sys.hcrit))
             )
             #Laplace pressure
             state.pressure[:,:,1] .-= (sys.gamma[2,3] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4]) .* 1/6 .* Swalbe.wetting_potential_93.(sys.hmin./(state.height[:,:,2] .+ sys.hcrit))) .* (
@@ -425,15 +425,15 @@ function filmpressure!(state::MultiLayer_2D, sys::SysConstMultiLayer)
 	elseif sys.n==3 && sys.m==2
       # if sys.layers == 2
       #Store what we will need twice
-      state.hi[:,:,1] .= (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_disj_32.(sys.hmin  ./(state.height[:,:,2] .+ sys.hcrit)) .* (1 .+ state.grad_h_sq ./2)
+      state.hi[:,:,1] .= (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_32.(sys.hmin  ./(state.height[:,:,2] .+ sys.hcrit)) .* (1 .+ state.grad_h_sq ./2)
         # disjoining pressure
         state.pressure[:,:,1] .=  (
-            (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])                  *(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_disj_32.(sys.hmin./(state.height[:,:,1] .+ sys.hcrit))
+            (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])                  *(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_32.(sys.hmin./(state.height[:,:,1] .+ sys.hcrit))
             .- state.hi[:,:,1]
         )
         state.pressure[:,:,2] .=  (
             state.hi[:,:,1]
-            .+ (sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_disj_32.(2*sys.hmin  ./(state.height[:,:,1] .+ state.height[:,:,2] + 2*sys.hcrit))
+            .+ (sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin)   .* Swalbe.fast_32.(2*sys.hmin  ./(state.height[:,:,1] .+ state.height[:,:,2] + 2*sys.hcrit))
         )
         #Laplace pressure
         state.pressure[:,:,1] .-= (
@@ -463,31 +463,6 @@ end
 
 
 
-function filmpressure!(state::StateActive, sys::SysConstActive)
-    hip, hjp, him, hjm, hipjp, himjp, himjm, hipjm = viewneighbors(state.dgrad)
-    # Straight elements j+1, i+1, i-1, j-1
-    circshift!(hip, state.height, (1,0))
-    circshift!(hjp, state.height, (0,1))
-    circshift!(him, state.height, (-1,0))
-    circshift!(hjm, state.height, (0,-1))
-    # Diagonal elements
-    circshift!(hipjp, state.height, (1,1))
-    circshift!(himjp, state.height, (-1,1))
-    circshift!(himjm, state.height, (-1,-1))
-    circshift!(hipjm, state.height, (1,-1))
-
-    # Stefan had a sign minus here, that isn't obvious to me. Them equation read positive
-    state.pressure .= -sys.γ_0 .* ((1 .- cospi(sys.θ_0)) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin)
-                      .* (Swalbe.power_broad.(sys.hmin./(state.height .+ sys.hcrit), sys.n)
-                      .- Swalbe.power_broad.(sys.hmin./(state.height .+ sys.hcrit), sys.m)) )
-    # Now the Laplace
-    state.pressure .-= sys.γ_0 .* (2/3 .* (hjp .+ hip .+ him .+ hjm)
-                   .+ 1/6 .* (hipjp .+ himjp .+ himjm .+ hipjm)
-                   .- 10/3 .* state.height)
-    return nothing
-end
-
-
 
 """
     surface_tension!(state::StateMiscible_1D, sys::SysConstMiscible_1D)
@@ -611,12 +586,12 @@ function filmpressure!(state::StateMiscible_1D, sys::SysConstMiscible_1D)
     circshift!(hip, state.height[:,1].+state.height[:,2], 1)
     circshift!(him, state.height[:,1].+state.height[:,2], -1)
 	if sys.n==9 && sys.m==3
-    		state.pressure[:,1] .= (-state.gamma[:,1] .- state.gamma[:,2] .+ state.gamma[:,3]) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_93.(sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ 2 .* sys.hcrit))
+    		state.pressure[:,1] .= (-state.gamma[:,1] .- state.gamma[:,2] .+ state.gamma[:,3]) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_93.(sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ 2 .* sys.hcrit))
             state.pressure[:,2] .=state.pressure[:,1]
             # For extra stability
 	    state.pressure .+= 0.1.* (-state.gamma[:,1] .- state.gamma[:,2] .+ state.gamma[:,3]) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin / 2) .* Swalbe.power_9.((sys.hmin/2)./(state.height .+ sys.hcrit))
 	elseif sys.n==3 && sys.m==2
-    		state.pressure .= (-state.gamma[:,1] .- state.gamma[:,2] .+ state.gamma[:,3]) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_32.(sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ sys.hcrit))
+    		state.pressure .= (-state.gamma[:,1] .- state.gamma[:,2] .+ state.gamma[:,3]) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_32.(sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ sys.hcrit))
 	else
         	throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
     	end
@@ -631,9 +606,9 @@ function filmpressure_fast!(state::Active_1D, sys::SysConstActive_1D)
     state.k .= ((-state.γ  .+ cospi.(sys.θ_0)*sys.γ_0 ) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin))
     #Π(h)
     if sys.n==9 && sys.m==3
-        state.pressure .= (state.k .* Swalbe.fast_disj_93.(sys.hmin ./ (state.height .+ sys.hcrit)))
+        state.pressure .= (state.k .* Swalbe.fast_93.(sys.hmin ./ (state.height .+ sys.hcrit)))
     elseif sys.n==3 && sys.m==2
-        state.pressure .= (state.k .* Swalbe.fast_disj_32.(sys.hmin ./ (state.height .+ sys.hcrit)))
+        state.pressure .= (state.k .* Swalbe.fast_32.(sys.hmin ./ (state.height .+ sys.hcrit)))
     else
         throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
     end
@@ -705,19 +680,19 @@ function filmpressure!(state::StateMultiLayer_1D, sys::SysConstMultiLayer_1D)
             state.pressure[:,1] .= - (hip[:,1] .- 2 .* state.height[:,1] .+ him[:,1])                                                .* (sys.gamma[2,3] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4]) .* 1/6 .* Swalbe.wetting_potential_93.(sys.hmin./(state.height[:,2] .+ sys.hcrit)))
             state.pressure[:,2] .= - (hip[:,1] .+ hip[:,2] .- 2 .* (state.height[:,1] .+ state.height[:,2]) .+ him[:,1] .+ him[:,2]) .*  sys.gamma[3,4]
             # store the one field we need twice
-            state.hi[:,1].=(sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*2*sys.hmin)   .* Swalbe.fast_disj_93.(2*sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit), sys.repulsive[1,2])
+            state.hi[:,1].=(sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*2*sys.hmin)   .* Swalbe.fast_93.(2*sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit))
             # disjoining pressure
-            state.pressure[:,1] .+= state.pressure[:,2] .+ state.hi[:,1] .+ (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_disj_93.(sys.hmin./(state.height[:,1] .+ sys.hcrit), sys.repulsive[1,1])
-            state.pressure[:,2] .+=                        state.hi[:,1] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_disj_93.(sys.hmin./(state.height[:,2] .+ sys.hcrit), sys.repulsive[2,1]).* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./2)
+            state.pressure[:,1] .+= state.pressure[:,2] .+ state.hi[:,1] .+ (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_93.(sys.hmin./(state.height[:,1] .+ sys.hcrit))
+            state.pressure[:,2] .+=                        state.hi[:,1] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_93.(sys.hmin./(state.height[:,2] .+ sys.hcrit)).* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./2)
         elseif sys.n==3 && sys.m==2
            #Laplace pressure
            state.pressure[:,1] .= - (hip[:,1] .- 2 .* state.height[:,1] .+ him[:,1])                                                .* (sys.gamma[2,3] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4]) .* Swalbe.wetting_potential_32.(sys.hmin./(state.height[:,2] .+ sys.hcrit)))
            state.pressure[:,2] .= - (hip[:,1] .+ hip[:,2] .- 2 .* (state.height[:,1] .+ state.height[:,2]) .+ him[:,1] .+ him[:,2]) .*  sys.gamma[3,4]
        # store the one field we need twice
-       state.hi[:,1].=(sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*2*sys.hmin)   .* Swalbe.fast_disj_32.(2*sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit), sys.repulsive[1,2])
+       state.hi[:,1].=(sys.gamma[2,3]+sys.gamma[1,4]-sys.gamma[2,4]-sys.gamma[1,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*2*sys.hmin)   .* Swalbe.fast_32.(2*sys.hmin./(state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit))
            # disjoining pressure
-           state.pressure[:,1] .+= state.pressure[:,2] .+ state.hi[:,1] .+ (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_disj_32.(  sys.hmin./(state.height[:,1] .+ sys.hcrit), sys.repulsive[1,1])
-           state.pressure[:,2] .+=                        state.hi[:,1] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_disj_32.(sys.hmin./(state.height[:,2] .+ sys.hcrit), sys.repulsive[2,1]).* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./2)
+           state.pressure[:,1] .+= state.pressure[:,2] .+ state.hi[:,1] .+ (sys.gamma[1,3]-sys.gamma[1,2]-sys.gamma[2,3])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_32.(  sys.hmin./(state.height[:,1] .+ sys.hcrit))
+           state.pressure[:,2] .+=                        state.hi[:,1] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4])*(sys.n-1)*(sys.m-1)/((sys.n-sys.m)*sys.hmin) .* Swalbe.fast_32.(sys.hmin./(state.height[:,2] .+ sys.hcrit)).* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./2)
         else
                 throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
         end
@@ -734,15 +709,15 @@ function filmpressure!(state::StateMultiLayer_1D, sys::SysConstMultiLayer_1D)
             # Store things we will need multiple times
             # When going to more then 3 layers one should really think about adding an extra field to store those especially because then we could use state.hi to store the actual interface positions
                 # \phi_123'(h_2)
-                state.hi[:,1] .= (sys.gamma[1,5] + sys.gamma[2,4] - sys.gamma[1,4] - sys.gamma[2,5]) * 1/3 * prefac * Swalbe.fast_disj_32.(3*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ state.height[:,3] .+ 3*sys.hcrit), sys.repulsive[1,3])
+                state.hi[:,1] .= (sys.gamma[1,5] + sys.gamma[2,4] - sys.gamma[1,4] - sys.gamma[2,5]) * 1/3 * prefac * Swalbe.fast_32.(3*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ state.height[:,3] .+ 3*sys.hcrit))
                 # \phi_{23}'(h_2+h_3)(1+(\nabla z_1)^2/2)
-                state.hi[:,2] .= (sys.gamma[2,5] + sys.gamma[3,4] - sys.gamma[2,4] -sys.gamma[3,5]) * 0.5 * prefac .* Swalbe.fast_disj_32.( 2*sys.hmin ./ (state.height[:,2] .+ state.height[:,3] .+ 2* sys.hcrit), sys.repulsive[2,2]) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
+                state.hi[:,2] .= (sys.gamma[2,5] + sys.gamma[3,4] - sys.gamma[2,4] -sys.gamma[3,5]) * 0.5 * prefac .* Swalbe.fast_32.( 2*sys.hmin ./ (state.height[:,2] .+ state.height[:,3] .+ 2* sys.hcrit)) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
                 # \phi_12'(h_3)(1+ (\nabla z_2)^2/2)
-                state.hi[:,3] .=(sys.gamma[1,4] + sys.gamma[2,3] - sys.gamma[1,3] - sys.gamma[2,4]) * 0.5 * prefac * Swalbe.fast_disj_32.(2*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit), sys.repulsive[1,2])
+                state.hi[:,3] .=(sys.gamma[1,4] + sys.gamma[2,3] - sys.gamma[1,3] - sys.gamma[2,4]) * 0.5 * prefac * Swalbe.fast_32.(2*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit))
             # Disjoining pressure
-                state.pressure[:,1] .+= state.pressure[:,2] .+ state.pressure[:,3] .+ state.hi[:,1]                  .+ state.hi[:,3] .+ (sys.gamma[1,3] - sys.gamma[1,2] - sys.gamma[2,3]) * prefac .* Swalbe.fast_disj_32.(sys.hmin ./ (state.height[:,1] .+ sys.hcrit), sys.repulsive[1,1])
-                state.pressure[:,2] .+= state.pressure[:,3]                        .+ state.hi[:,1] .+ state.hi[:,2] .+ state.hi[:,3] .+ (sys.gamma[2,4] - sys.gamma[2,3] - sys.gamma[3,4]) * prefac .* Swalbe.fast_disj_32.(sys.hmin ./ (state.height[:,2] .+ sys.hcrit), sys.repulsive[2,1]) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
-                state.pressure[:,3] .+=                                               state.hi[:,1] .+ state.hi[:,2]                  .+ (sys.gamma[3,5] - sys.gamma[3,4] - sys.gamma[4,5]) * prefac .* Swalbe.fast_disj_32.(sys.hmin ./ (state.height[:,3] .+ sys.hcrit), sys.repulsive[3,1]) .* (1 .+ state.grad_h[:,2] .* state.grad_h[:,2] ./ 2)
+                state.pressure[:,1] .+= state.pressure[:,2] .+ state.pressure[:,3] .+ state.hi[:,1]                  .+ state.hi[:,3] .+ (sys.gamma[1,3] - sys.gamma[1,2] - sys.gamma[2,3]) * prefac .* Swalbe.fast_32.(sys.hmin ./ (state.height[:,1] .+ sys.hcrit))
+                state.pressure[:,2] .+= state.pressure[:,3]                        .+ state.hi[:,1] .+ state.hi[:,2] .+ state.hi[:,3] .+ (sys.gamma[2,4] - sys.gamma[2,3] - sys.gamma[3,4]) * prefac .* Swalbe.fast_32.(sys.hmin ./ (state.height[:,2] .+ sys.hcrit)) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
+                state.pressure[:,3] .+=                                               state.hi[:,1] .+ state.hi[:,2]                  .+ (sys.gamma[3,5] - sys.gamma[3,4] - sys.gamma[4,5]) * prefac .* Swalbe.fast_32.(sys.hmin ./ (state.height[:,3] .+ sys.hcrit)) .* (1 .+ state.grad_h[:,2] .* state.grad_h[:,2] ./ 2)
             elseif sys.n==9 && sys.m==3
                  # Laplace pressure
                  state.pressure[:,1] .= - (hip[:,1] .- 2 .* state.height[:,1] .+ him[:,1])                                                                                             .* (sys.gamma[2,3] .+ (sys.gamma[2,4]-sys.gamma[2,3]-sys.gamma[3,4]) .* 1/6 .* Swalbe.wetting_potential_93.(sys.hmin./(state.height[:,2] .+ sys.hcrit)) .+ (sys.gamma[1,5] + sys.gamma[3,4] - sys.gamma[2,4] - sys.gamma[3,5]) * 1/6 .* Swalbe.wetting_potential_93.(2*sys.hmin ./ (state.height[:,2] .+ state.height[:,3] .+ 2*sys.hcrit)))
@@ -753,14 +728,14 @@ function filmpressure!(state::StateMultiLayer_1D, sys::SysConstMultiLayer_1D)
              # Store things we will need multiple times
              # When going to more then 3 layers one should really think about adding an extra field to store those especially because then we could use state.hi to store the actual interface positions
                  # \phi_123'(h_2)
-                 state.hi[:,1] .= (sys.gamma[1,5] + sys.gamma[2,4] - sys.gamma[1,4] - sys.gamma[2,5]) * 1/3 * prefac *fast_disj_93.(3*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ state.height[:,3] .+ 3*sys.hcrit), sys.repulsive[1,3])
+                 state.hi[:,1] .= (sys.gamma[1,5] + sys.gamma[2,4] - sys.gamma[1,4] - sys.gamma[2,5]) * 1/3 * prefac *Swalbe.fast_93.(3*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ state.height[:,3] .+ 3*sys.hcrit))
                  # \phi_{23}'(h_2+h_3)(1+(\nabla z_1)^2/2)
-                 state.hi[:,2] .= (sys.gamma[2,5] + sys.gamma[3,4] - sys.gamma[2,4] -sys.gamma[3,5]) * 0.5 * prefac .*fast_disj_93.( 2*sys.hmin ./ (state.height[:,2] .+ state.height[:,3] .+ 2* sys.hcrit), sys.repulsive[2,2]) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
+                 state.hi[:,2] .= (sys.gamma[2,5] + sys.gamma[3,4] - sys.gamma[2,4] -sys.gamma[3,5]) * 0.5 * prefac .*Swalbe.fast_93.( 2*sys.hmin ./ (state.height[:,2] .+ state.height[:,3] .+ 2* sys.hcrit)) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
                  # \phi_12'(h_3)(1+ (\nabla z_2)^2/2)
-                 state.hi[:,3] .=(sys.gamma[1,4] + sys.gamma[2,3] - sys.gamma[1,3] - sys.gamma[2,4]) * 0.5 * prefac *fast_disj_93.(2*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit), sys.repulsive[1,2])
+                 state.hi[:,3] .=(sys.gamma[1,4] + sys.gamma[2,3] - sys.gamma[1,3] - sys.gamma[2,4]) * 0.5 * prefac *Swalbe.fast_93.(2*sys.hmin ./ (state.height[:,1] .+ state.height[:,2] .+ 2*sys.hcrit))
              # Disjoining pressure
-                 state.pressure[:,1] .+= state.pressure[:,2] .+ state.pressure[:,3] .+ state.hi[:,1]                  .+ state.hi[:,3] .+ (sys.gamma[1,3] - sys.gamma[1,2] - sys.gamma[2,3]) * prefac .*fast_disj_93.(sys.hmin ./ (state.height[:,1] .+ sys.hcrit), sys.repulsive[1,1])
-                 state.pressure[:,2] .+= state.pressure[:,3]                        .+ state.hi[:,1] .+ state.hi[:,2] .+ state.hi[:,3] .+ (sys.gamma[2,4] - sys.gamma[2,3] - sys.gamma[3,4]) * prefac .*fast_disj_93.(sys.hmin ./ (state.height[:,2] .+ sys.hcrit), sys.repulsive[2,1]) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
+                 state.pressure[:,1] .+= state.pressure[:,2] .+ state.pressure[:,3] .+ state.hi[:,1]                  .+ state.hi[:,3] .+ (sys.gamma[1,3] - sys.gamma[1,2] - sys.gamma[2,3]) * prefac .*Swalbe.fast_93.(sys.hmin ./ (state.height[:,1] .+ sys.hcrit))
+                 state.pressure[:,2] .+= state.pressure[:,3]                        .+ state.hi[:,1] .+ state.hi[:,2] .+ state.hi[:,3] .+ (sys.gamma[2,4] - sys.gamma[2,3] - sys.gamma[3,4]) * prefac .*Swalbe.fast_93.(sys.hmin ./ (state.height[:,2] .+ sys.hcrit)) .* (1 .+ state.grad_h[:,1] .* state.grad_h[:,1] ./ 2)
         else
                 throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
         end
@@ -829,15 +804,15 @@ function filmpressure_curved!(state::State_curved_1D, sys::SysConst_1D)
     circshift!(hip, state.height .+ state.substrate, 1)
     circshift!(him, state.height .+ state.substrate, -1)
         if sys.n==9 && sys.m==3
-                state.pressure .= -(1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate) .* sys.γ .* (1 .- cospi.(sys.θ)) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_93.(sys.hmin./(state.height .+ sys.hcrit))
-        elseif sys.n==3 && sys.m==2
-                state.pressure .= - (1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate) .* sys.γ .* (1 .- cospi.(sys.θ)) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_32.(sys.hmin./(state.height .+ sys.hcrit))
-        elseif sys.n==0 && sys.m==0
+                state.pressure .= -(1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate) .* sys.param.γ .* (1 .- cospi.(sys.param.θ)) .* (sys.param.n - 1) .* (sys.param.m - 1) ./ ((sys.param.n - sys.param.m) * sys.param.hmin) .* Swalbe.fast_93.(sys.param.hmin./(state.height .+ sys.param.hcrit))
+        elseif sys.param.n==3 && sys.param.m==2
+                state.pressure .= - (1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate) .* sys.param.γ .* (1 .- cospi.(sys.param.θ)) .* (sys.param.n - 1) .* (sys.param.m - 1) ./ ((sys.param.n - sys.param.m) * sys.param.hmin) .* Swalbe.fast_32.(sys.param.hmin./(state.height .+ sys.param.hcrit))
+        elseif sys.param.n==0 && sys.param.m==0
                 state.pressure .= 0.0
         else
-                throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
+                throw(DomainError((sys.param.n,sys.param.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
         end
-        state.pressure .-= sys.γ .* (hip .- 2 .* (state.height .+ state.substrate) .+ him)
+        state.pressure .-= sys.param.γ .* (hip .- 2 .* (state.height .+ state.substrate) .+ him)
     return nothing
 end
 
@@ -881,15 +856,15 @@ function filmpressure_curved_healing!(state::State_curved_1D, sys::SysConst_1D,S
 	# Straight elements j+1, i+1, i-1, j-1
     circshift!(hip, state.height .+ state.substrate, 1)
     circshift!(him, state.height .+ state.substrate, -1)
-        if sys.n==9 && sys.m==3
-                state.pressure .= S .* (1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate)* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_93.(sys.hmin./(state.height .+ sys.hcrit))
-        elseif sys.n==3 && sys.m==2                                                                                                                                         state.pressure .= S .* (1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_32.(sys.hmin./(state.height .+ sys.hcrit))
-        elseif sys.n==0 && sys.m==0
+        if sys.param.n==9 && sys.param.m==3
+                state.pressure .= S .* (1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate)* (sys.param.n - 1) .* (sys.param.m - 1) ./ ((sys.param.n - sys.param.m) * sys.param.hmin) .* Swalbe.fast_93.(sys.param.hmin./(state.height .+ sys.param.hcrit))
+        elseif sys.param.n==3 && sys.param.m==2                                                                                                                                         state.pressure .= S .* (1 .+ 0.5 .* state.grad_substrate .* state.grad_substrate) .* (sys.param.n - 1) .* (sys.param.m - 1) ./ ((sys.param.n - sys.param.m) * sys.param.hmin) .* Swalbe.fast_32.(sys.param.hmin./(state.height .+ sys.param.hcrit))
+        elseif sys.param.n==0 && sys.param.m==0
                 state.pressure .= 0.0
         else
-		throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
+		throw(DomainError((sys.param.n,sys.param.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
 	end
-	state.pressure .-= sys.γ .* (hip .- 2 .* (state.height .+ state.substrate) .+ him)
+	state.pressure .-= sys.param.γ .* (hip .- 2 .* (state.height .+ state.substrate) .+ him)
     return nothing
 end
 
diff --git a/test/fancy_models.jl b/test/fancy_models.jl
index dd3826a..99bd279 100644
--- a/test/fancy_models.jl
+++ b/test/fancy_models.jl
@@ -9,7 +9,7 @@
     		tdump=Int(floor(Tmax/10)), # dump time
     		h= 1,                       # initial film height
     		rho=1,                      # initial catalyst concentration
-    		ϵ=0.001,                    # initial noise
+    		eps=0.001,                    # initial noise
     		γ_0=0.001,                  # reference surface tension in the absence of chemical compounds
     		Γ = 0.1,                    # Surface tension effect product rho_A
     		GammaB=0.1,                 # Surface tension effect reactant rhoB
@@ -61,7 +61,7 @@
         		state.rho_A .= rho_0A
         		state.rho_B .= rho_0B
         		# This set up the film height as constant plus a small noise without high wavenumber frequencies
-        		Swalbe.browniannoise!(state.height, h, ϵ, L/2-L/10)
+			state.height .= h .* ones(L) .* eps .* rand(L)
 			height_0 = zeros(L)
 			height_0 .= state.height
         			# The LBM time loop
@@ -95,7 +95,7 @@
 	@testset "Curved Substrates" begin 
 	println("testing liquid films on curved substrates")
 		function curved(;
-        		Tmax=1e3,                       # Simulation time
+        		Tmax=1000,                       # Simulation time
         		tdump=Int(floor(Tmax/10)),     # dump time
         		h= 1,                           # initial film height
         		eps=0.001,                      # initial noise
@@ -107,7 +107,7 @@
         		omega=3,                        # number of periods of the underlying substrate
 		)
         		# set up system
-        		sys = Swalbe.SysConst_1D{Float64}(L =  L ,Tmax=Tmax, tdump=tdump, γ = gamma, δ=delta,θ =theta, n=n,m=m,hmin=hmin)
+			sys=Swalbe.SysConst_1D(L=L, param=Swalbe.Taumucs(Tmax=Tmax, tdump=tdump, γ = gamma, δ=delta,θ =theta, n=n,m=m,hmin=hmin))
         		state = Swalbe.Sys(sys, kind="curved")
         		# Initial data
         		state.height.= h .+ eps .* rand(sys.L)
diff --git a/test/runtests.jl b/test/runtests.jl
index b624fc6..a72e7de 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,26 +1,26 @@
 using Swalbe, Test, Statistics, Parameters, Random, LazySets, BSON, FileIO, DataFrames
 
 @testset "Swalbe.jl" begin
-    println("Testing allocations")
-    include("initialize.jl")
-    println("Testing initial states")
-    include("initialvalues.jl")
-    println("Testing finite differences")
-    include("differences.jl")
-    println("Testing film pressure")
-    include("pressure.jl")
-    println("Testing collision operator")
-    include("collide.jl")
-    println("Testing equilibrium distribution")
-    include("equilibrium.jl")
-    println("Testing forces")
-    include("forcing.jl")
-    println("Testing moment calculation")
-    include("moments.jl")
-    println("Testing measurements")
-    include("measures.jl")
-    println("Testing with relevant simulations")
-    include("simulate.jl")
+    # println("Testing allocations")
+    # include("initialize.jl")
+    # println("Testing initial states")
+    # include("initialvalues.jl")
+    # println("Testing finite differences")
+    # include("differences.jl")
+    # println("Testing film pressure")
+    # include("pressure.jl")
+    # println("Testing collision operator")
+    # include("collide.jl")
+    # println("Testing equilibrium distribution")
+    # include("equilibrium.jl")
+    # println("Testing forces")
+    # include("forcing.jl")
+    # println("Testing moment calculation")
+    # include("moments.jl")
+    # println("Testing measurements")
+    # include("measures.jl")
+    # println("Testing with relevant simulations")
+    # include("simulate.jl")
     println("Testing models beyond flat films")
     include("fancy_models.jl")
 end