diff --git a/fpp-ptc-sandbox/code/z_radiation_matt_fake_maps.f90 b/fpp-ptc-sandbox/code/z_radiation_matt_fake_maps.f90
index 32406c2..e53251d 100644
--- a/fpp-ptc-sandbox/code/z_radiation_matt_fake_maps.f90
+++ b/fpp-ptc-sandbox/code/z_radiation_matt_fake_maps.f90
@@ -135,13 +135,13 @@ program example
    
   do i=1,c_%nd2
   m%v(i)=m%v(i)*decrement(i)
-  enddo
- m1=m + m
+ enddo
+ !m1=3.d0 + m
 
-   do i=1,6
-    write(6,format6) real(m1%e_ij(i,1:6) )
-  enddo
-  stop
+  ! do i=1,6
+  !  write(6,format6) real(m1%e_ij(i,1:6) )
+  !enddo
+  !stop
 
 ! do i=1,1000
 !  m1  = m1*m1
diff --git a/src/NonlinearNormalForm.jl b/src/NonlinearNormalForm.jl
index 40d9149..7d4ad4a 100644
--- a/src/NonlinearNormalForm.jl
+++ b/src/NonlinearNormalForm.jl
@@ -19,7 +19,9 @@ import Base: ∘,
              copy!,
              convert,
              show,
-             rand
+             rand,
+             promote_rule,
+             eltype
 
 import LinearAlgebra: norm,
                       dot,
diff --git a/src/map/compose.jl b/src/map/compose.jl
index 2e0bfcd..abe4f61 100644
--- a/src/map/compose.jl
+++ b/src/map/compose.jl
@@ -28,8 +28,7 @@ See the documentation for `compose_it!` for information on `work_low` and `work_
 - `work_prom` -- Temporary vector of allocated `ComplexTPS`s when there is implicit promotion. See the `compose_it!` documentation for more details. Default is output from `prep_comp_work_prom(m, m2, m1)`
 """
 function compose!(m::DAMap, m2::DAMap, m1::DAMap; keep_scalar::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, dospin::Bool=true, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_comp_work_low(m), work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
-  checkop(m, m2, m1)
-  checkpromotion(m, m2, m1)
+  checkinplace(m, m2, m1)
   
   # DAMap setup:
   desc = getdesc(m1)
@@ -92,8 +91,7 @@ See the documentation for `compose_it!` for information on `work_low` and `work_
 - `work_prom` -- Temporary vector of allocated `ComplexTPS`s when there is implicit promotion. See the `compose_it!` documentation for more details. Default is output from `prep_comp_work_prom(m, m2, m1)`
 """
 function compose!(m::TPSAMap, m2::TPSAMap, m1::TPSAMap; dospin::Bool=true, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_comp_work_low(m), work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
-  checkop(m, m2, m1)
-  checkpromotion(m, m2, m1)
+  checkinplace(m, m2, m1)
   
   # TPSAMap setup:
   # For TPSA Map concatenation, we need to subtract w_0 (m2 x0) (Eq. 33)
@@ -134,41 +132,7 @@ $($t) composition, which calculates `m2 ∘ m1` $( $t == DAMap ? "ignoring the s
 """
 function compose(m2::$t,m1::$t)
   checkop(m2, m1)
-
-  desc = getdesc(m1)
-  nn = numnn(desc)
-  nv = numvars(desc)
-
-  outT = promote_type(eltype(m2.x),eltype(m1.x))
-  
-  # set up outx0
-  outx0 = Vector{numtype(outT)}(undef, nv)
-
-  # Set up outx:
-  outx = Vector{outT}(undef, nn)
-  for i=1:nv  # no need to allocate immutable parameters taken care of inside compose_it!
-      @inbounds outx[i] = outT(use=desc)
-  end
-
-  # set up quaternion out:
-  if !isnothing(m1.Q)
-    outq = Vector{outT}(undef, 4)
-    for i=1:4
-      @inbounds outq[i] = outT(use=desc)
-    end
-    outQ = Quaternion(outq)
-  else
-    outQ = nothing
-  end
-
-  # set up stochastic out
-  if isnothing(m1.E) && isnothing(m2.E)
-    outE = nothing
-  else
-    outE = Matrix{numtype(outT)}(undef, nv, nv)
-  end
-
-  m = $t(outx0, outx, outQ, outE, m1.idpt)
+  m = zero_op(m1,m2)
   compose!(m, m2, m1)
   
   return m
diff --git a/src/map/compose_it.jl b/src/map/compose_it.jl
index 488b380..71b0df5 100644
--- a/src/map/compose_it.jl
+++ b/src/map/compose_it.jl
@@ -42,9 +42,7 @@ Note that the `ComplexTPS`s in the vector(s) must be allocated and have the same
 If spin is included, not that the final quaternion concatenation step mul! will creat allocations
 """ 
 function compose_it!(m::$t, m2::$t, m1::$t; dospin::Bool=true, dostochastic::Bool=true, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_comp_work_low(m), work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
-  checkop(m, m2, m1)
-  checkpromotion(m, m2, m1)
- 
+  checkinplace(m, m2, m1)
   @assert !(m === m1) "Cannot compose_it!(m, m2, m1) with m === m1"
 
   desc = getdesc(m1)
@@ -153,16 +151,8 @@ function compose_it!(m::$t, m2::$t, m1::$t; dospin::Bool=true, dostochastic::Boo
   # Stochastic
   # MAKE THIS FASTER!
   if !isnothing(m.E) && dostochastic
-    if !isnothing(m1.E)
-      M2 = jacobian(m2)   
-      m.E .= M2*m1.E*transpose(M2) 
-    else 
-      m.E .= 0
-    end
-
-    if !isnothing(m2.E)
-      m.E .+= m2.E
-    end
+    M2 = jacobian(m2)   
+    m.E .= M2*m1.E*transpose(M2) + m2.E
   end
 
   return 
diff --git a/src/map/ctors.jl b/src/map/ctors.jl
index bce11b5..d89cfef 100644
--- a/src/map/ctors.jl
+++ b/src/map/ctors.jl
@@ -2,7 +2,7 @@ for t = (:DAMap, :TPSAMap)
 @eval begin
 
 """
-    $($t)(m::TaylorMap{S,T,U,V,W}; use::UseType=nothing, idpt::Union{Nothing,Bool}=m.idpt) where {S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
+    $($t)(m::Union{TaylorMap{S,T,U,V,W},Probe{S,T,U,V,W}}; use=nothing, idpt::Union{Nothing,Bool}=m.idpt) where {S,T,U,V,W}
 
 Creates a new copy of the passed `TaylorMap` as a `$($t)`. 
 
@@ -13,152 +13,54 @@ parameters must agree, however the orders may be different.
 
 If `idpt` is not specified, then the same `idpt` as `m` will be used, else that specified will be used.
 """
-function $t(m::TaylorMap{S,T,U,V,W}; use::UseType=nothing, idpt::Union{Nothing,Bool}=m.idpt) where {S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
-  if isnothing(use)
-    use = getdesc(m)
-  else
-    numnn(use) == numnn(m) || error("Number of variables + parameters in GTPSAs for `m` and `use` disagree!")
-  end
+function $t(m::Union{TaylorMap{S,T,U,V,W},Probe{S,T,U,V,W}}; use=m, idpt::Union{Nothing,Bool}=m.idpt) where {S,T<:Union{TPS,ComplexTPS},U,V,W}
+  numnn(use) == numnn(m) || error("Number of variables + parameters in GTPSAs for `m` and `use` disagree!")
   
+  outm = $t{S,T,U,V,typeof(idpt)}(undef, use = use, idpt = idpt) # undefined map but parameters allocated
   nv = numvars(use)
-  np = numparams(use)
-  nn = numnn(use)
 
-
-  x = Vector{T}(undef, nn)
+  # allocate variables
   for i=1:nv
-    x[i] = T(m.x[i], use=getdesc(use))
+    @inbounds outm.x[i] = T(m.x[i], use=getdesc(use))
   end
 
-  # use same parameters if same descriptor (use=nothing)
-  if isnothing(use) || getdesc(use) == getdesc(m)
-    @inbounds x[nv+1:nn] .= view(m.x, nv+1:nn)
-  else
-    if T == TPS
-      @inbounds x[nv+1:nn] .= params(getdesc(first(x)))
-    else
-      @inbounds x[nv+1:nn] .= complexparams(getdesc(first(x)))
-    end
-  end
-
-  if !isnothing(m.Q)
-    q = Vector{T}(undef, 4)
+  # allocate quaternion if present
+  if !isnothing(outm.Q)
     for i=1:4
-      @inbounds q[i] = T(m.Q.q[i],use=getdesc(use))
+      @inbounds outm.Q.q[i] = T(m.Q.q[i],use=getdesc(use))
     end
-    Q = Quaternion(q)
-  else
-    Q = nothing
   end
 
-  x0 = Vector{S}(undef, nv)
+  # set the reference orbit properly
   if nv > numvars(m)  # Increasing dimensionality
-    x0[1:numvars(m)] .= m.x0
-    x0[numvars(m)+1:nv] .= 0
+    outm.x0[1:numvars(m)] .= m.x0
+    outm.x0[numvars(m)+1:nv] .= 0
   else    # Reducing or keeping same dimensionality
-    x0[1:nv] .= view(m.x0, 1:nv)
+    outm.x0[1:nv] .= view(m.x0, 1:nv)
   end
 
-  if !isnothing(m.E)
-    E = similar(m.E, nv, nv)
+  # set the FD matrix properly
+  if !isnothing(outm.E)
     if nv > numvars(m)  # Increasing dimensionality
-      E[1:numvars(m),1:numvars(m)] .= view(m.E, 1:numvars(m), 1:numvars(m))
-      E[numvars(m)+1:nv,:] .= 0
-      E[:,numvars(m)+1:nv] .= 0
-    else
-      E[1:nv,1:nv] .= view(m.E, 1:nv, 1:nv)
-    end
-  else
-    E = nothing
-  end
-
-  return $t{S,T,U,V,typeof(idpt)}(x0, x, Q, E, idpt)
-end
-
-"""
-    $($t)(p::Probe{S,T,U,V,W}; use::UseType=nothing, idpt::Union{Nothing,Bool}=p.idpt) where {S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
-
-Creates a `$($t)` from the `Probe`, which must contain `TPS`s. 
-
-If `use` is not specified, then the same GTPSA `Descriptor` as `p` will be used. If `use` is 
-specified (could be another `Descriptor`, `TaylorMap`, or a `Probe` containing `TPS`s), then the 
-`p` promoted to a $($(t)) will have the same `Descriptor` as in `use.` The total number of variables + 
-and parameters must agree, however the orders may be different.
-
-If `idpt` is not specified, then the same `idpt` as `m` will be used, else that specified will be used.
-"""
-function $t(p::Probe{S,T,U,V,W}; use::UseType=nothing, idpt::Union{Nothing,Bool}=p.idpt) where {S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
-  if isnothing(use)
-    use = getdesc(p)
-  else
-    numnn(use) == numnn(p) || error("Number of variables + parameters in GTPSAs for `p` and `use` disagree!")
-  end
-  
-  nv = numvars(use)
-  np = numparams(use)
-  nn = numnn(use)
-  
-  x = Vector{T}(undef, nn)
-
-  length(p.x) == numvars(p) || error("Length of orbital ray ($(length(p.x))) inconsistent with number of variables in GTPSA ($(nv))")
-  
-  for i=1:numvars(p)
-    @inbounds x[i] = T(p.x[i], use=getdesc(use))
-  end
-
-  for i=numvars(p)+1:nv
-    @inbounds x[i] = T(use=getdesc(use))
-  end
-
-  if T == TPS
-    @inbounds x[nv+1:nn] .= params(getdesc(first(x)))
-  else
-    @inbounds x[nv+1:nn] .= complexparams(getdesc(first(x)))
-  end
-
-
-  if !isnothing(p.Q)
-    q = Vector{T}(undef, 4)
-    for i=1:4
-      @inbounds q[i] = T(p.Q.q[i],use=getdesc(use))
-    end
-    Q = Quaternion(q)
-  else
-    Q = nothing
-  end
-
-  x0 = Vector{S}(undef, nv)
-  if nv > numvars(p)  # Increasing dimensionality
-    x0[1:numvars(p)] .= p.x0
-    x0[numvars(p)+1:nv] .= 0
-  else    # Reducing or keeping same dimensionality
-    x0[1:nv] .= view(p.x0, 1:nv)
-  end
-
-  if !isnothing(p.E)
-    E = similar(p.E, nv, nv)
-    if nv > numvars(p)  # Increasing dimensionality
-      E[1:numvars(p),1:numvars(p)] .= view(p.E, 1:numvars(p), 1:numvars(p))
-      E[numvars(p)+1:nv,:] .= 0
-      E[:,numvars(p)+1:nv] .= 0
+      outm.E[1:numvars(m),1:numvars(m)] .= view(m.E, 1:numvars(m), 1:numvars(m))
+      outm.E[numvars(m)+1:nv,:] .= 0
+      outm.E[:,numvars(m)+1:nv] .= 0
     else
-      E[1:nv,1:nv] .= view(p.E, 1:nv, 1:nv)
+      outm.E[1:nv,1:nv] .= view(m.E, 1:nv, 1:nv)
     end
-  else
-    E = nothing
   end
 
-  return $t{S,T,U,V,typeof(idpt)}(x0, x, Q, E, idpt)
+  return outm
 end
 
 """
-    $($t){S,T,U,V,W}(u::UndefInitializer; use::UseType=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
+    $($t){S,T,U,V,W}(u::UndefInitializer; use=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
 
 Creates an undefined `$($t){S,T,U,V,W}` with same `Descriptor` as `use`. The immutable 
 parameters will be allocated if `use` is not a `TaylorMap`, else the immutable parameters 
 from `use` will be used.
 """
-function $t{S,T,U,V,W}(u::UndefInitializer; use::UseType=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
+function $t{S,T,U,V,W}(u::UndefInitializer; use=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
   desc = getdesc(use)
   nv = numvars(desc)
   np = numparams(desc)
@@ -166,6 +68,8 @@ function $t{S,T,U,V,W}(u::UndefInitializer; use::UseType=GTPSA.desc_current, idp
   x0 = Vector{S}(undef, nv)
   x = Vector{T}(undef, nn)
 
+  
+
   # use same parameters if use isa TaylorMap and T == eltype(use.x)
   if use isa TaylorMap && T == eltype(use.x)
     @inbounds x[nv+1:nn] .= view(use.x, nv+1:nn)
@@ -178,130 +82,128 @@ function $t{S,T,U,V,W}(u::UndefInitializer; use::UseType=GTPSA.desc_current, idp
   end
 
   if U != Nothing
-    q = Vector{T}(undef, 4)
+    q = Vector{eltype(U)}(undef, 4)
     Q = Quaternion(q)
   else
     Q = nothing
   end
 
   if V != Nothing
-    E = Matrix{S}(undef, nv, nv)
+    E = Matrix{eltype(V)}(undef, nv, nv)
   else
     E = nothing
   end
 
-  return $t(x0, x, Q, E, idpt)
+  return $t{S,T,U,V,W}(x0, x, Q, E, idpt)
 end
 
 """
-    $($t)(; x::Union{Vector{<:Union{TPS,ComplexTPS}},Nothing}=nothing, x0::Union{Vector,Nothing}=nothing, Q::Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing}=nothing, E::Union{Matrix,Nothing}=nothing,  spin::Union{Bool,Nothing}=nothing, stochastic::Union{Bool,Nothing}=nothing, idpt::Union{Nothing,Bool}=nothing, use::UseType=nothing)
+    $($t)(; x::Union{Vector{<:Union{TPS,ComplexTPS}},Nothing}=nothing, x0::Union{Vector,Nothing}=nothing, Q::Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing}=nothing, E::Union{Matrix,Nothing}=nothing,  spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing, idpt::Union{Nothing,Bool}=nothing, use=nothing)
 
 Constructs a $($t) with the passed vector of `TPS`/`ComplexTPS` as the orbital ray, and optionally the entrance 
-coordinates `x0`, `Quaternion` for spin `Q`, and stochastic matrix `E` as keyword arguments. The helper keyword 
-arguments `spin` and `stochastic` may be set to `true` to construct a $($t) with an identity quaternion/stochastic 
-matrix, or `false` for no spin/stochastic. Note that setting `spin`/`stochastic` to any `Bool` value without `Q` or `E` 
+coordinates `x0`, `Quaternion` for spin `Q`, and FD matrix `E` as keyword arguments. The helper keyword 
+arguments `spin` and `FD` may be set to `true` to construct a $($t) with an identity quaternion/FD 
+matrix, or `false` for no spin/FD. Note that setting `spin`/`FD` to any `Bool` value without `Q` or `E` 
 specified is type-unstable. This constructor also checks for consistency in the length of the orbital ray and GTPSA 
 `Descriptor`. The `use` kwarg may also be used to change the `Descriptor` of the TPSs, provided the number of variables 
 + parameters agree (orders may be different).
 """
-function $t(; x::Union{Vector{<:Union{TPS,ComplexTPS}},Nothing}=nothing, x0::Union{Vector,Nothing}=nothing, Q::Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing}=nothing, E::Union{Matrix,Nothing}=nothing,  spin::Union{Bool,Nothing}=nothing, stochastic::Union{Bool,Nothing}=nothing, idpt::Union{Nothing,Bool}=nothing, use::UseType=nothing)
-  if isnothing(use)
-    use = GTPSA.desc_current
+function $t(;use=GTPSA.desc_current, x::Vector{T}=vars(getdesc(use)), x0::Vector{S}=zeros(numtype(eltype(x)), numvars(use)), Q::U=nothing, E::V=nothing, idpt::W=nothing, spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing) where {S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
+  # set up
+  if isnothing(spin)
+    if isnothing(Q)
+      outU = Nothing
+    else
+      outU = typeof(Q)
+    end
+  elseif spin
+    if isnothing(Q)
+      outU = Quaternion{T}
+    else
+      outU = typeof(Q)
+    end
+  else
+    error("For no spin tracking, please omit the spin kwarg or set spin=nothing") # For type stability
+    #outU = Nothing # For type instability
+  end
+  
+  if isnothing(FD)
+    if isnothing(E)
+      outV = Nothing
+    else
+      outV = typeof(E)
+    end
+  elseif FD
+    if isnothing(E)
+      outV = Matrix{numtype(T)}
+    else
+      outV = typeof(E)
+    end
+  else
+    error("For no fluctuation-dissipation, please omit the FD kwarg or set FD=nothing") # For type stability
+    #outV = Nothing # For type instability
   end
 
+  outm = $t{S,T,outU,outV,typeof(idpt)}(undef, use = use, idpt = idpt)   
+
   nv = numvars(use)
   np = numparams(use)
   nn = numnn(use)
 
-  if isnothing(x)
-    x1 = Vector{TPS}(undef, nn)
-    for i=1:nv
-      @inbounds x1[i] = TPS(use=use)
-    end
-  else
-    numvars(use) == numvars(x) || error("Number of variables in GTPSAs for `x` and `use` disagree!")
-    numvars(x) == length(x) || error("Length of orbital ray inconsistent with number of variables in GTPSA!")
-    x1 = Vector{eltype(x)}(undef, nn)
-    @inbounds x1[1:nv] .= view(x, 1:nv)
-  end
+  # sanity checks
+  length(x0) == nv || error("Number of variables != length of reference orbit vector!")
+  length(x) == nv || error("Number of variables in GTPSAs for `x` and `use` disagree!")
 
-  if eltype(x1) == TPS
-    @inbounds x1[nv+1:nn] .= params(getdesc(use))
-  else
-    @inbounds x1[nv+1:nn] .= complexparams(getdesc(use))
-  end
 
-  if isnothing(x0)
-    x01 = zeros(numtype(first(x1)), nv)
-  else
-    numvars(x1) == length(x0) || error("Number of variables != length of reference orbit vector!")
-    x01 = copy(x0)
+  outm.x0 .= x0
+
+  for i=1:nv
+    @inbounds outm.x[i] = T(x[i], use=getdesc(use))
   end
 
-  if isnothing(spin)
-    if isnothing(Q)
-      Q1 = Q
-    else
-      (!isnothing(use) || getdesc(x) == getdesc(Q)) || error("Orbital ray Descriptor different from quaternion Descriptor!")
-      q = Vector{eltype(x1)}(undef, 4)
+  if !isnothing(outm.Q)
+    if !isnothing(Q)
       for i=1:4
-        @inbounds q[i] = (eltype(x1))(Q.q[i],use=getdesc(use))
+        @inbounds outm.Q.q[i] = eltype(outU)(Q.q[i], use=getdesc(use))
       end
-      Q1 = Quaternion(q)
-    end
-  elseif spin
-    if isnothing(Q)
-      Q1 = Quaternion(first(x1)) # implicilty uses use descriptor
     else
-      (!isnothing(use) || getdesc(x) == getdesc(Q)) || error("Orbital ray Descriptor different from quaternion Descriptor!")
-      q = Vector{eltype(x1)}(undef, 4)
       for i=1:4
-        @inbounds q[i] = (eltype(x1))(Q.q[i],use=getdesc(use))
+        @inbounds outm.Q.q[i] = (eltype(outU))(use=getdesc(use))
       end
-      Q1 = Quaternion(q)
+      outm.Q.q[1][0] = 1
     end
-  else
-    error("For no spin tracking, please omit the spin kwarg or set spin=nothing") # For type stability
-    #Q1 = nothing # For type instability
   end
 
-  if isnothing(stochastic)
-    E1 = E
-  elseif stochastic
-    if isnothing(E)
-      E1 = zeros(eltype(x01), nv, nv) 
+  if !isnothing(outm.E)
+    if !isnothing(E)
+      (nv,nv) == size(E) || error("Size of FD matrix inconsistent with number of variables!")
+      outm.E .= E
     else
-      (nv,nv) == size(E) || error("Size of stochastic matrix inconsistent with number of variables!")
-      E1 = E
+      outm.E .= 0
     end
-  else
-    error("For no stochastic, please omit the stochastic kwarg or set stochastic=nothing") # For type stability
-    #E1 = nothing # for type instability
   end
 
-  return $t(x01, x1, Q1, E1, idpt)
+  return outm
 end
 
 """
-    $($t)(M; use::UseType=GTPSA.desc_current, x0::Vector{S}=zeros(eltype(M), size(M,1)), Q::U=nothing, E::V=nothing,  spin::Union{Bool,Nothing}=nothing, stochastic::Union{Bool,Nothing}=nothing, idpt::Union{Nothing,Bool}=nothing) where {S,U<:Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing},V<:Union{Matrix,Nothing}}
+    $($t)(M::AbstractMatrix; use=GTPSA.desc_current, x0::Vector{S}=zeros(numtype(T), numvars(use)), Q::U=nothing, E::V=nothing, idpt::W=nothing, spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing) where {S,U<:Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
 
 `M` must represent a matrix with linear indexing.
 
 Constructs a $($t) with the passed matrix of scalars `M` as the linear part of the `TaylorMap`, and optionally the entrance 
-coordinates `x0`, `Quaternion` for spin `Q`, and stochastic matrix `E` as keyword arguments. The helper keyword 
-arguments `spin` and `stochastic` may be set to `true` to construct a $($t) with an identity quaternion/stochastic 
-matrix, or `false` for no spin/stochastic. Note that setting `spin`/`stochastic` to any `Bool` value without `Q` or `E` 
+coordinates `x0`, `Quaternion` for spin `Q`, and FD matrix `E` as keyword arguments. The helper keyword 
+arguments `spin` and `FD` may be set to `true` to construct a $($t) with an identity quaternion/FD 
+matrix, or `false` for no spin/FD. Note that setting `spin`/`FD` to any `Bool` value without `Q` or `E` 
 specified is type-unstable. This constructor also checks for consistency in the length of the orbital ray and GTPSA 
 `Descriptor`.
 """
-function $t(M; use::UseType=GTPSA.desc_current, x0::Vector{S}=zeros(eltype(M), numvars(use)), Q::U=nothing, E::V=nothing,  spin::Union{Bool,Nothing}=nothing, stochastic::Union{Bool,Nothing}=nothing, idpt::Union{Nothing,Bool}=nothing) where {S,U<:Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing},V<:Union{Matrix,Nothing}}
+function $t(M::AbstractMatrix; use=GTPSA.desc_current, x::Vector{T}=vars(getdesc(use)), x0::Vector{S}=zeros(numtype(eltype(x)), numvars(use)), Q::U=nothing, E::V=nothing, idpt::W=nothing, spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing) where {S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{<:Union{TPS,ComplexTPS}},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
   Base.require_one_based_indexing(M)
   nv = numvars(use)
-  np = numparams(use)
   nn = numnn(use)
 
-  nv == length(x0) || error("Number of variables in GTPSA Descriptor != length of reference orbit vector!")
-  nv >= size(M,1) || error("Number of rows in transfer matrix > number of variables in GTPSA!")
+  nv >= size(M,1) || error("Number of rows in matrix > number of variables in GTPSA!")
 
   if eltype(M) <: Complex
     outT = ComplexTPS
@@ -309,107 +211,30 @@ function $t(M; use::UseType=GTPSA.desc_current, x0::Vector{S}=zeros(eltype(M), n
     outT = TPS
   end
 
-  x1 = Vector{outT}(undef, nn)
+  x = Vector{outT}(undef, nv)
   for i=1:size(M,1)
-    @inbounds x1[i] = (outT)(use=getdesc(use))
+    @inbounds x[i] = (outT)(use=getdesc(use))
     for j=1:size(M,2)
-      @inbounds x1[i][j] = M[i,j]
+      @inbounds x[i][j] = M[i,j]
     end
   end
 
-  for i=size(M,1)+1:nv
-    @inbounds x1[i] = (outT)(use=getdesc(use))
-  end
-
-  if outT == TPS
-    @inbounds x1[nv+1:nn] .= params(getdesc(first(x1)))
-  else
-    @inbounds x1[nv+1:nn] .= complexparams(getdesc(first(x1)))
-  end
-
-  if isnothing(spin)
-    if isnothing(Q)
-      Q1 = Q
-    else
-      q = Vector{outT}(undef, 4)
-      for i=1:4
-        @inbounds q[i] = outT(Q.q[i],use=getdesc(use))
-      end
-      Q1 = Quaternion(q)
-    end
-  elseif spin
-    if isnothing(Q)
-      Q1 = Quaternion(first(x1)) # implicilty uses use descriptor
-    else
-      q = Vector{outT}(undef, 4)
-      for i=1:4
-        @inbounds q[i] = outT(Q.q[i],use=getdesc(use))
-      end
-      Q1 = Quaternion(q)
-    end
-  else
-    error("For no spin tracking, please omit the spin kwarg or set spin=nothing") # For type stability
-    #Q1 = nothing # For type instability
-  end
-
-  if isnothing(stochastic)
-    E1 = E
-  elseif stochastic
-    if isnothing(E)
-      E1 = zeros(eltype(x01), nv, nv) 
-    else
-      (nv,nv) == size(E) || error("Size of stochastic matrix inconsistent with number of variables!")
-      E1 = E
-    end
-  else
-    error("For no stochastic, please omit the stochastic kwarg or set stochastic=nothing") # For type stability
-    #E1 = nothing # for type instability
-  end
-  return $t{eltype(M),outT,typeof(Q1),typeof(E1),typeof(idpt)}(copy(x0), x1, Q1, E1,idpt)
+  return DAMap(use=use, x=x, x0=x0,Q=Q,E=E,idpt=idpt,spin=spin,FD=FD)
 end
 
 
 """
     zero(m::$($t))
 
-Creates a $($t) with the same GTPSA `Descriptor`, and spin/stochastic on/off,
+Creates a $($t) with the same GTPSA `Descriptor`, and spin/FD on/off,
 as `m` but with all zeros for each quantity (except for the immutable parameters 
 in `x[nv+1:nn]`, which will be copied from `m.x`)
 """
 function zero(m::$t)
-  desc = getdesc(m)
-  nn = numnn(desc)
-  nv = numvars(desc)
-  np = numparams(desc)
-  
-  x = Vector{eltype(m.x)}(undef, nn)
-  for i=1:nv
-    @inbounds x[i] = (eltype(m.x))(use=desc)
-  end
-
-  # use same parameters 
-  @inbounds x[nv+1:nn] .= view(m.x, nv+1:nn)
-
-  if !isnothing(m.Q)
-    q = Vector{eltype(m.x)}(undef, 4)
-    for i=1:4
-      @inbounds q[i] = (eltype(m.x))(use=desc)
-    end
-    Q = Quaternion(q)
-  else
-    Q = nothing
-  end
-
-  if !isnothing(m.E)
-    E = zeros(eltype(m.E), nv, nv)
-  else
-    E = nothing
-  end
-
-  return $t(zeros(eltype(m.x0), nv), x, Q, E, m.idpt)
+  return zero(typeof(m), use=m, idpt=m.idpt)
 end
 
-function zero(::Type{$t{S,T,U,V,W}}; use::UseType=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
+function zero(::Type{$t{S,T,U,V,W}}; use=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
   desc = getdesc(use)
   nn = numnn(desc)
   nv = numvars(desc)
@@ -435,9 +260,9 @@ function zero(::Type{$t{S,T,U,V,W}}; use::UseType=GTPSA.desc_current, idpt::W=no
   end
 
   if U != Nothing
-    q = Vector{T}(undef, 4)
+    q = Vector{eltype(U)}(undef, 4)
     for i=1:4
-      @inbounds q[i] = T(use=desc)
+      @inbounds q[i] = eltype(U)(use=desc)
     end
     Q = Quaternion(q)
   else
@@ -445,109 +270,95 @@ function zero(::Type{$t{S,T,U,V,W}}; use::UseType=GTPSA.desc_current, idpt::W=no
   end
 
   if V != Nothing
-    E = zeros(S, nv, nv)
+    E = Matrix{eltype(V)}(undef, nv, nv) #zeros(eltype(V), nv, nv)
   else
     E = nothing
   end
 
-  return $t(x0, x, Q, E, idpt)
+  return $t{S,T,U,V,W}(x0, x, Q, E, idpt)
+end
+
+function zero_op(m2::Union{$t,Number}, m1::Union{$t,Number})
+  outtype = promote_type(typeof(m1),typeof(m2))
+
+  # If either inputs are ComplexTPS, use those parameters
+  if m2 isa $t
+    if m1 isa $t
+      if eltype(m2.x) == ComplexTPS
+        return zero(outtype,use=m2,idpt=m2.idpt)
+      else
+        return zero(outtype,use=m1,idpt=m1.idpt)
+      end
+    else
+      return zero(outtype,use=m2,idpt=m2.idpt)
+    end
+  elseif m1 isa $t
+    return zero(outtype,use=m1,idpt=m1.idpt)
+  else
+    error("Cannot create new map based only on scalars")
+  end
 end
 
+
 """
     one(m::$($t))
   
 Construct an identity map based on `m`.
 """
 function one(m::$t)
-  desc = getdesc(m)
-  nn = numnn(desc)
-  nv = numvars(desc)
-  np = numparams(desc)
+  return one(typeof(m), use=m, idpt=m.idpt)
+end
 
-  T = eltype(m.x)
+function one(t::Type{$t{S,T,U,V,W}}; use=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
+  m = zero(t, use=use, idpt=idpt)
+  nv = numvars(m)
   
-  x = Vector{T}(undef, nn)
-  if T == ComplexTPS
-    for i=1:nv
-      @inbounds x[i] = complexmono(i,use=desc)
-    end
-  else
-    for i=1:nv
-      @inbounds x[i] = mono(i,use=desc)
-    end
+  for i=1:nv
+    @inbounds m.x[i][i] = 1
   end
 
-  # use same parameters 
-  @inbounds x[nv+1:nn] .= view(m.x, nv+1:nn)
-
   if !isnothing(m.Q)
-    q = Vector{T}(undef, 4)
-    @inbounds q[1] = one(first(x))
-    for i=2:4
-      @inbounds q[i] = T(use=desc)
-    end
-    Q = Quaternion(q)
-  else
-    Q = nothing
+    @inbounds m.Q.q[1][0] = 1
   end
 
-  if !isnothing(m.E)
-    E = zeros(eltype(m.E), nv, nv)
-  else
-    E = nothing
-  end
-
-  return $t(zeros(eltype(m.x0), nv), x, Q, E, m.idpt)
+  return m
 end
 
-function one(::Type{$t{S,T,U,V,W}}; use::UseType=GTPSA.desc_current, idpt::W=nothing) where {S,T,U,V,W}
-  desc = getdesc(use)
+end
+end
+#=
+function test1(m2,m1)
+  desc = getdesc(m1)
   nn = numnn(desc)
   nv = numvars(desc)
-  np = numparams(desc)
 
-  x0 = zeros(S, nv)
+  outT = promote_type(eltype(m2.x),eltype(m1.x))
 
-  x = Vector{T}(undef, nn)
-  for i=1:nv
-    @inbounds x[i] = T(use=desc)
-    x[i][i] = 1
-  end
+  # set up outx0
+  outx0 = Vector{numtype(outT)}(undef, nv)
 
-  if use isa Union{TaylorMap,Probe} && eltype(use.x) == T
-    # use same parameters 
-    @inbounds x[nv+1:nn] .= view(m.x, nv+1:nn)
-  else
-    # allocate
-    if T == TPS
-      @inbounds x[nv+1:nn] .= params(desc)
-    else
-      @inbounds x[nv+1:nn] .= complexparams(desc)
-    end
+  # Set up outx:
+  outx = Vector{outT}(undef, nn)
+  for i=1:nv  # no need to allocate immutable parameters taken care of inside compose_it!
+      @inbounds outx[i] = outT(use=desc)
   end
 
-  if U != Nothing
-    q = Vector{T}(undef, 4)
+  # set up quaternion out:
+  if !isnothing(m1.Q)
+    outq = Vector{outT}(undef, 4)
     for i=1:4
-      @inbounds q[i] = T(use=desc)
+      @inbounds outq[i] = outT(use=desc)
     end
-    q[1][0] = 1
-    Q = Quaternion(q)
+    outQ = Quaternion(outq)
   else
-    Q = nothing
+    outQ = nothing
   end
 
-  if V != Nothing
-    E = zeros(S, nv, nv)
+  # set up FD out
+  if isnothing(m1.E) && isnothing(m2.E)
+    outE = nothing
   else
-    E = nothing
+    outE = Matrix{numtype(outT)}(undef, nv, nv)
   end
-
-  return $t(x0, x, Q, E, idpt)
-end
-
-end
-end
-
-
-
+  return DAMap(outx0, outx, outQ, outE, m1.idpt)
+end=#
\ No newline at end of file
diff --git a/src/map/inv.jl b/src/map/inv.jl
index ecc78bf..45dfa37 100644
--- a/src/map/inv.jl
+++ b/src/map/inv.jl
@@ -39,7 +39,7 @@ end
 =#
 
 """
-    inv!(m::TaylorMap{S,T,U,V,W}, m1::TaylorMap{S,T,U,V,W}; dospin::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_inv_work_low(m1)) where {S,T,U,V}
+    inv!(m::TaylorMap, m1::TaylorMap; dospin::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_inv_work_low(m1))
 
 In-place inversion of the `TaylorMap` setting `m = inv(m1)`. Aliasing `m === m1` is allowed, however 
 in this case a temporary vector must be used to store the scalar part of `m1` prior to inversion so 
@@ -50,8 +50,8 @@ that the entrance/exit coordinates of the map can be properly handled.
 - `work_ref` -- If `m === m1`, then a temporary vector must be used to store the scalar part. If not provided and `m === m1`, this temporary will be created internally. Default is `nothing`
 - `work_low` -- Temporary vector to hold the low-level C pointers. Default is output from `prep_inv_work_low`
 """
-function inv!(m::TaylorMap{S,T,U,V,W}, m1::TaylorMap{S,T,U,V,W}; dospin::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_inv_work_low(m1)) where {S,T,U,V,W}
-  checkop(m, m1)
+function inv!(m::TaylorMap, m1::TaylorMap; dospin::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_inv_work_low(m1))
+  checkinplace(m, m1)
 
   desc = getdesc(m1)
   nn = numnn(desc)
diff --git a/src/map/methods.jl b/src/map/methods.jl
index 59581c2..b0f5339 100644
--- a/src/map/methods.jl
+++ b/src/map/methods.jl
@@ -68,8 +68,7 @@ end
 
 # --- copy! ---
 function copy!(m::TaylorMap, m1::TaylorMap)
-  checkidpt(m, m1)
-  checkspin(m, m1)
+  checkinplace(m, m1)
 
   m.x0 .= m1.x0
   desc = getdesc(m)
@@ -122,16 +121,15 @@ jacobiant(m::TaylorMap;include_params=false) = GTPSA.jacobiant(view(m.x, 1:numva
 checksymp(m::TaylorMap) = checksymp(GTPSA.jacobian(m))
 
 # --- cutord ---
-function cutord(m1::TaylorMap{S,T,U,V}, order::Integer, spin_order::Integer=order; dospin::Bool=true) where {S,T,U,V}
+function cutord(m1::TaylorMap, order::Integer, spin_order::Integer=order; dospin::Bool=true)
   m = zero(m1)
   cutord!(m, m1, order, spin_order, dospin=dospin)
   return m
 end
 
-function cutord!(m::TaylorMap{S,T,U,V}, m1::TaylorMap{S,T,U,V}, order::Integer, spin_order::Integer=order; dospin::Bool=true) where {S,T,U,V}
-  checkidpt(m, m1)
-  checkspin(m, m1)
-  
+function cutord!(m::TaylorMap, m1::TaylorMap, order::Integer, spin_order::Integer=order; dospin::Bool=true)
+  checkinplace(m, m1)
+
   desc = getdesc(m1)
   nv = numvars(desc)
   np = numparams(desc)
@@ -155,15 +153,14 @@ function cutord!(m::TaylorMap{S,T,U,V}, m1::TaylorMap{S,T,U,V}, order::Integer,
 end
 
 # --- getord ---
-function getord(m1::TaylorMap{S,T,U,V}, order::Integer, spin_order::Integer=order; dospin::Bool=true) where {S,T,U,V}
+function getord(m1::TaylorMap, order::Integer, spin_order::Integer=order; dospin::Bool=true)
   m = zero(m1)
   getord!(m, m1, order, spin_order, dospin=dospin)
   return m
 end
 
-function getord!(m::TaylorMap{S,T,U,V}, m1::TaylorMap{S,T,U,V}, order::Integer, spin_order::Integer=order; dospin::Bool=true) where {S,T,U,V}
-  checkidpt(m, m1)
-  checkspin(m, m1)
+function getord!(m::TaylorMap, m1::TaylorMap, order::Integer, spin_order::Integer=order; dospin::Bool=true)
+  checkinplace(m, m1)
   
   desc = getdesc(m1)
   nv = numvars(desc)
@@ -185,4 +182,20 @@ function getord!(m::TaylorMap{S,T,U,V}, m1::TaylorMap{S,T,U,V}, order::Integer,
     m.E .= m.E
   end
   return
+end
+
+# --- setmatrix! ---
+
+function setmatrix!(m::DAMap, M::AbstractMatrix)
+  Base.require_one_based_indexing(M)
+  nv = numvars(m)
+  nn = numnn(m)
+
+  nv >= size(M,1) || error("Number of rows in matrix > number of variables in GTPSA!")
+
+  for i=1:size(M,1)
+    for j=1:size(M,2)
+      @inbounds m.x[i][j] = M[i,j]
+    end
+  end
 end
\ No newline at end of file
diff --git a/src/map/operators.jl b/src/map/operators.jl
index c61c859..fb78661 100644
--- a/src/map/operators.jl
+++ b/src/map/operators.jl
@@ -115,10 +115,7 @@ for ops = (("add!", :+), ("sub!",:-))
 @eval begin
 
 function $(Meta.parse(ops[1]))(m::TaylorMap, m1::TaylorMap, m2::TaylorMap; dospin::Bool=true)
-  checkidpt(m, m1, m2)
-  checkspin(m, m1, m2)
-
-  assertstoch(m, m1, m2)
+  checkinplace(m, m1, m2)
   
   nv = numvars(m)
 
@@ -148,10 +145,7 @@ function $(Meta.parse(ops[1]))(m::TaylorMap, m1::TaylorMap, m2::TaylorMap; dospi
 end
 
 function $(Meta.parse(ops[1]))(m::TaylorMap, J::UniformScaling, m1::TaylorMap; dospin::Bool=true)
-  checkidpt(m, m1)
-  checkspin(m, m1)
-  
-  assertstoch(m,m1,m)
+  checkinplace(m, m1)
   
   nv = numvars(m)
 
@@ -181,10 +175,7 @@ function $(Meta.parse(ops[1]))(m::TaylorMap, J::UniformScaling, m1::TaylorMap; d
 end
 
 function $(Meta.parse(ops[1]))(m::TaylorMap, m1::TaylorMap, J::UniformScaling; dospin::Bool=true)
-  checkidpt(m, m1)
-  checkspin(m, m1)
-
-  assertstoch(m,m1,m)
+  checkinplace(m, m1)
 
   nv = numvars(m)
 
@@ -212,15 +203,9 @@ function $(Meta.parse(ops[1]))(m::TaylorMap, m1::TaylorMap, J::UniformScaling; d
 end
 
 function $(ops[2])(m1::TaylorMap, m2::TaylorMap)
-  checkidpt(m1, m2)
-  checkspin(m1, m2)
+  checkop(m1, m2)
+  m = zero_op(m1, m2)
 
-  # Promote if necessary:
-  if eltype(m1.x) == ComplexTPS
-    m = zero(m1)
-  else
-    m = zero(m2)
-  end
   $(Meta.parse(ops[1]))(m, m1, m2)
   return m
 end
@@ -246,8 +231,7 @@ end
 for ops = (("add!", :+), ("sub!",:-), ("mul!",:*), ("div!",:/))
 @eval begin
 function $(Meta.parse(ops[1]))(m::TaylorMap, a::Number, m1::TaylorMap; dospin::Bool=true)
-  checkidpt(m, m1)
-  checkspin(m, m1)
+  checkinplace(m, a, m1)
   
   nv = numvars(m)
   m.x0 .= m1.x0
@@ -269,8 +253,7 @@ function $(Meta.parse(ops[1]))(m::TaylorMap, a::Number, m1::TaylorMap; dospin::B
 end
 
 function $(Meta.parse(ops[1]))(m::TaylorMap, m1::TaylorMap, a::Number; dospin::Bool=true)
-  checkidpt(m, m1)
-  checkspin(m, m1)
+  checkinplace(m, a, m1)
   
   nv = numvars(m)
 
@@ -293,21 +276,15 @@ function $(Meta.parse(ops[1]))(m::TaylorMap, m1::TaylorMap, a::Number; dospin::B
 end
 
 function $(ops[2])(a::Number, m1::TaylorMap)
-  if a isa Complex
-    m = zero(complex(typeof(m1)), use=m1)
-  else
-    m = zero(m1)
-  end
+  m = zero_op(m1, a)
+
   $(Meta.parse(ops[1]))(m, a, m1)
   return m
 end
 
 function $(ops[2])(m1::TaylorMap, a::Number)
-  if a isa Complex
-    m = zero(complex(typeof(m1)), use=m1)
-  else
-    m = zero(m1)
-  end
+  m = zero_op(m1, a)
+  
   $(Meta.parse(ops[1]))(m, m1, a)
   return m
 end
diff --git a/src/normal.jl b/src/normal.jl
index e95dcbc..18295f8 100644
--- a/src/normal.jl
+++ b/src/normal.jl
@@ -36,7 +36,9 @@ function normal(m::DAMap)
     a1_inv_matrix[nhv+1,nhv+1] = 1; a1_inv_matrix[nhv+2,nhv+2] = 1;
   end
 
-  a1 = DAMap(complex(inv(a1_inv_matrix)),use=m0,idpt=m.idpt,E=zeros(ComplexF64, numvars(m),numvars(m)))
+  a1 = zero(promote_type(eltype(a1_inv_matrix),typeof(m0)),use=m0,idpt=m.idpt)
+  setmatrix!(a1, inv(a1_inv_matrix))
+
   m1 = a1^-1 ∘ m0 ∘ a1
 
   # 3: Go into phasor's basis
@@ -120,7 +122,7 @@ end
 
 function equilibrium_moments(m::DAMap, a::DAMap)
   !isnothing(m.E) || error("Map does not have stochasticity")
-  !all(m.E .== 0) || error("No stochastic fluctuations in map (m.E .== 0)")
+  !all(m.E .== 0) || error("No FD fluctuations in map (m.E .== 0)")
 
   # Moments Σ transform with map like MΣMᵀ + B 
   # This is only linear because this is very complex to higher order not necessary
@@ -128,9 +130,9 @@ function equilibrium_moments(m::DAMap, a::DAMap)
   # To include parameters later, we would:
   # First must go to parameter-dependent fixed point to all orders (obtained from factorization after 
   # normal form) and then can get the a1 matrix around there
-  # tracking code would have to give stochastic matrix as a function of the parameters
+  # tracking code would have to give FD matrix as a function of the parameters
 
-  # Let B = m.E (stochastic part)
+  # Let B = m.E (FD part)
   # We want to find Σ such that Σ = MΣMᵀ + B  (fixed point)
   # very easy to do in phasors basis
   # fixed point transformation does nothing (note a0.E = Bₐ₀ = 0 of course and a0.x is identity in variable but not in parameters)
diff --git a/src/probe.jl b/src/probe.jl
index dff5eb4..6466638 100644
--- a/src/probe.jl
+++ b/src/probe.jl
@@ -1,4 +1,4 @@
-function Probe(x::Vector{T}; x0::Vector{S}=zeros(length(x)), Q::U=nothing, E::V=nothing, spin::Union{Bool,Nothing}=nothing, stochastic::Union{Bool,Nothing}=nothing, idpt::W=nothing) where {S,T,U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
+function Probe(x::Vector{T}; x0::Vector{S}=zeros(length(x)), Q::U=nothing, E::V=nothing, spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing, idpt::W=nothing) where {S,T,U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
   length(x) == length(x0) || error("Length of orbital ray != length of reference orbit vector!")
 
   if isnothing(spin)
@@ -14,17 +14,17 @@ function Probe(x::Vector{T}; x0::Vector{S}=zeros(length(x)), Q::U=nothing, E::V=
     #Q1 = nothing # For type instability
   end
 
-  if isnothing(stochastic)
+  if isnothing(FD)
     E1 = E
-  elseif stochastic
+  elseif FD
     if isnothing(E)
       E1 = zeros(eltype(x0), length(x), length(x)) 
     else
-      (length(x),length(x)) == size(E) || error("Size of stochastic matrix inconsistent with number of variables!")
+      (length(x),length(x)) == size(E) || error("Size of FD matrix inconsistent with number of variables!")
       E1 = E
     end
   else
-    error("For no stochasticity, please omit the stochastic kwarg or set stochastic=nothing") # For type stability
+    error("For no stochasticity, please omit the FD kwarg or set FD=nothing") # For type stability
     #E1 = nothing # for type instability
   end
 
diff --git a/src/types.jl b/src/types.jl
index b3a564c..61f4358 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,13 +1,13 @@
 
 """
-    TaylorMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}}
+    TaylorMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
 
 Abstract type for `TPSAMap` and `DAMap` used for normal form analysis. 
 
 All `TaylorMap`s contain `x0` and `x` as the entrance coordinates and transfer map 
 as a truncated power series respectively. If spin is included, a field `Q` containing 
 a `Quaternion` as a truncated power series is included, else `Q` is `nothing`. If 
-stochasticity is included, a field `E` contains a matrix of the envelope for the stochastic
+stochasticity is included, a field `E` contains a matrix of the envelope for the FD
 kicks, else `E` is nothing.
 
 If all planes are exhibiting pseudo-harmonic oscillations, then `idpt` is `nothing`. 
@@ -19,18 +19,18 @@ is constant (energy-like): `idpt=false` if the variable with index `NV-1` is con
 - `x0`   -- Entrance coordinates of the map, Taylor expansion point
 - `x`    -- Orbital ray as a truncated power series, expansion around `x0` + scalar part equal to EXIT coordinates of map
 - `Q`    -- `Quaternion` as a truncated power series if spin is included, else `nothing`
-- `E`    -- Matrix of the envelope for stochastic kicks if included, else `nothing`
+- `E`    -- Matrix of the envelope for FD kicks if included, else `nothing`
 - `idpt` -- If the last plane is coasting, then `idpt=false` if the first variable in last plane is constant (energy-like) or `true` if the second. If no coasting, then `idpt=nothing`
 """
-abstract type TaylorMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}} end 
+abstract type TaylorMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}} end 
 
 """
-    DAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
+    DAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
 
 `TaylorMap` that composes and inverses as a `DAMap` (with the scalar part ignored).
 See `TaylorMap` for more information.
 """
-struct DAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
+struct DAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
   x0::Vector{S}    # Entrance value of map
   x::Vector{T}     # Expansion around x0, with scalar part equal to EXIT value of map wrt initial coordinates x0
   Q::U             # Quaternion for spin
@@ -39,12 +39,12 @@ struct DAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union
 end
 
 """
-    TPSAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
+    TPSAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
 
 `TaylorMap` that composes and inverses as a `TPSAMap` (with the scalar part included).
 See `TaylorMap` for more information.
 """
-struct TPSAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
+struct TPSAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}} <: TaylorMap{S,T,U,V,W}
   x0::Vector{S}    # Entrance value of map
   x::Vector{T}     # Expansion around x0, with scalar part equal to EXIT value of map wrt initial coordinates x0
   Q::U             # Quaternion for spin
@@ -53,19 +53,19 @@ struct TPSAMap{S,T<:Union{TPS,ComplexTPS},U<:Union{Quaternion{T},Nothing},V<:Uni
 end
 
 """
-    Probe{S,T,U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}}
+    Probe{S,T,U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
 
 Parametric type used for tracking. The orbital/spin part can contain 
 either scalars or `TPS`s. If spin is included, the field `Q` is a `Quaternion` 
 with the same type as `x`, else `Q` is `nothing`. If radiation is included, 
-the field `E` contains the stochastic matrix with `eltype(E)==S`, else `E` is `nothing`.
+the field `E` contains the FD matrix with `eltype(E)==S`, else `E` is `nothing`.
 
 If all planes are exhibiting pseudo-harmonic oscillations, then `idpt` is `nothing`. 
 If the last plane is coasting, then `idpt` specifies which variable in the last plane 
 is constant (energy-like): `idpt=false` if the variable with index `NV-1` is constant, or 
 `idpt=true` is the variable with index `NV` is constant.
 """
-struct Probe{S,T,U<:Union{Quaternion{T},Nothing},V<:Union{Matrix{S},Nothing},W<:Union{Nothing,Bool}}
+struct Probe{S,T,U<:Union{Quaternion{T},Nothing},V<:Union{Matrix,Nothing},W<:Union{Nothing,Bool}}
   x0::Vector{S}   # Entrance coordinates
   x::Vector{T}    # Out coordinates
   Q::U            # Quaternion
@@ -85,11 +85,29 @@ end
 
 
 
-UseType = Union{Descriptor, TPS, ComplexTPS, TaylorMap, Probe{<:Any,Union{TPS,ComplexTPS},<:Any,<:Any}, VectorField, Nothing}
-#=
-function promote_rule(m1::Type{TaylorMap{S1,T1,U,V1,W}}, m2::Type{TaylorMap{S2,T2,U,V2,W}}) where {S1,S2,T1,T2,U,V1,V2,W}
-  S = promote_rule(eltype(m1.x0), eltype(m2.x0))
-  T = promote_rule(eltype(m1.x), eltype(m2.x))
-  if !isnothing
-end =#
+#UseType = Union{Descriptor, TPS, ComplexTPS, DAMap, TPSAMap, Probe{<:Any,Union{TPS,ComplexTPS},<:Any,<:Any}, VectorField, Nothing}
+
+for t = (:DAMap, :TPSAMap)
+@eval begin    
+
+function promote_rule(::Type{$t{S,T,U,V,W}}, ::Type{G}) where {S,T,U,V,W,G<:Union{Number,Complex}}
+  outS = promote_type(S,numtype(T),G)
+  outT = promote_type(T,G)
+  U != Nothing ? outU = Quaternion{promote_type(eltype(U), G)} : outU = Nothing
+  V != Nothing ? outV = promote_type(Matrix{G},V) : outV = Nothing
+  return $t{outS,outT,outU,outV,W}
+end
+
+function promote_rule(::Type{$t{S1,T1,U1,V1,W}}, ::Type{$t{S2,T2,U2,V2,W}}) where {S1,S2,T1,T2,U1,U2,V1,V2,W} 
+  outS = promote_type(S1, S2, numtype(T1), numtype(T2))
+  outT = promote_type(T1, T2)
+  U1 != Nothing ? outU = Quaternion{promote_type(eltype(U2),eltype(U2))} : outU = Nothing
+  V1 != Nothing ? outV = promote_type(V1,V2) : outV = Nothing
+  return $t{outS,outT,outU,outV,W}
+end 
+
+end
+end
+
+
 
diff --git a/src/utils/misc.jl b/src/utils/misc.jl
index 7ef5c9d..a1f3b6c 100644
--- a/src/utils/misc.jl
+++ b/src/utils/misc.jl
@@ -24,42 +24,55 @@ vpords(m::Union{Probe{<:Real,<:Union{TPS,ComplexTPS},<:Any,<:Any},<:TaylorMap,Ve
 vords(m::Union{Probe{<:Real,<:Union{TPS,ComplexTPS},<:Any,<:Any},<:TaylorMap,VectorField}) = unsafe_wrap(Vector{UInt8}, unsafe_load(getdesc(m).desc).no, numvars(m))
 pords(m::Union{Probe{<:Real,<:Union{TPS,ComplexTPS},<:Any,<:Any},<:TaylorMap,VectorField}) = unsafe_wrap(Vector{UInt8}, unsafe_load(getdesc(m).desc).no, numparams(m))
 
-@inline checkidpt(m::TaylorMap...) = all(x->x.idpt==first(m).idpt, m) || error("Maps have disagreeing idpt")
-@inline checkspin(m...) = all(x->isnothing(x.Q), m) || all(x->!isnothing(x.Q), m) || error("Atleast one map/vector field includes spin while others do not")
+@inline checkidpt(maps::TaylorMap...) = all(x->x.idpt==first(maps).idpt, maps) || error("Maps have disagreeing idpt")
+@inline checkspin(stuff...) = all(x->isnothing(x.Q), stuff) || all(x->!isnothing(x.Q), stuff) || error("Atleast one map/vector field includes spin while others do not")
+@inline checkFD(maps::TaylorMap...) = all(x->isnothing(x.E), maps) || all(x->!isnothing(x.E), maps) || error("Atleast one map includes fluctuations/dissipation while others do not")
 
-@inline checkop(m::TaylorMap...) = checkidpt(m...) && checkspin(m...)
+@inline checkop(stuff...) = checkidpt(filter(x->(x isa TaylorMap), stuff)...) && checkspin(filter(x->(x isa Union{TaylorMap,VectorField}), stuff)...) && checkFD(filter(x->(x isa TaylorMap), stuff)...)
+
+# checkinplace is preferred to using the "where {S,..}.." syntax as this also ensure idpt is consistent
+# and checks promotion as well. It also gives descriptive errors rather than just "function not found"
+@inline function checkinplace(m::Union{TaylorMap,VectorField}, stuff...)
+  checkop(m, stuff...)
 
-@inline function checkdestination(m::TaylorMap, maps::TaylorMap...)
   # Checks that the output map has all types properly promoted
-  mapx0types = map(x->eltype(x.x0), maps)
-  outx0type = promote_type(mapx0types...)
-  eltype(m.x0) == outx0type || error("Output map reference orbit type $(eltype(m.x0)) must be $outx0type")
+  # or NOT promoted if unneccessary
+  maps = filter(x->(x isa TaylorMap), stuff)
+  mapsvfs = filter(x->(x isa Union{TaylorMap,VectorField}), stuff)
+  nums = filter(x->(x isa Number), stuff)
+  numtypes = map(x->typeof(x), nums) # scalars only affect x and Q, not x0 or E in FPP
+
+  xtypes = map(x->eltype(x.x), mapsvfs)
+  xnumtypes = map(x->numtype(x),xtypes)
+
+  if m isa TaylorMap
+    x0types = map(x->eltype(x.x0), maps)
+    outx0type = promote_type(x0types..., xnumtypes...) # reference orbit in composition is affected by orbital part
+    eltype(m.x0) == outx0type || error("Output $(typeof(m)) reference orbit type $(eltype(m.x0)) must be $outx0type")
+  end
 
-  mapxtypes = map(x->numtype(eltype(x.x)),maps)
-  outxtype = promote_type(mapxtypes...)
-  eltype(m.x) == outxtype || error("Output map orbital ray type $(eltype(m.x)) must be $outxtype")
+  outxtype = promote_type(xtypes..., numtypes...)
+  eltype(m.x) == outxtype || error("Output $(typeof(m)) orbital ray type $(eltype(m.x)) must be $xtype")
 
   if !isnothing(m.Q)
-    mapQtypes = map(x->numtype(eltype(x.Q.q)), maps)
-    outQtype = promote_type(mapQtypes...)
-    eltype(m.Q) == outQtype || error("Output map quaternion type $(eltype(m.Q.q)) must be $outQtype")
+    qtypes = map(x->numtype(eltype(x.Q.q)), mapsvfs)
+    outqtype = promote_type(qtypes..., numtypes)
+    eltype(m.Q) == outqtype || error("Output $(typeof(m)) quaternion type $(eltype(m.Q.q)) must be $outqtype")
   end
 
   # Part of the promotion is stochasticity:
   # the output map must include stochasticity if any input includes stochasticity:
-  !isnothing(m.E) || all(x->isnothing(x.E), maps) || error("Output map must include stochastic matrix (at least 1 input map includes stochastic matrix)")
-
-  if isnothing(m.E)
-    maptypes = map(x->numtype(eltype(x.x)),maps)
-    stochtypes = map(x->isnothing(x.E) ? Float64 : eltype(x.E), maps)
-    outtype = promote_type(maptypes..., stochtypes...)
-    eltype(m.E) == outtype || error("Output map stochastic matrix type $(eltype(m.E)) must be $outtype")
+  if m isa TaylorMap && !isnothing(m.E)
+    Etypes = map(x->eltype(x.E), maps)
+    outtype = promote_type(xnumtypes..., Etypes...)
+    eltype(m.E) == outtype || error("Output $(typeof(m)) FD matrix type $(eltype(m.E)) must be $outtype")
   end
+
   return true
 end
 
 # --- random symplectic map ---
-function rand(t::Union{Type{DAMap},Type{TPSAMap}}; spin::Union{Bool,Nothing}=nothing, stochastic::Union{Bool,Nothing}=nothing, use::Union{Descriptor,TPS,ComplexTPS}=GTPSA.desc_current, ndpt::Union{Nothing,Integer}=nothing)
+function rand(t::Union{Type{DAMap},Type{TPSAMap}}; spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing, use::Union{Descriptor,TPS,ComplexTPS}=GTPSA.desc_current, ndpt::Union{Nothing,Integer}=nothing)
   if isnothing(spin)
     U = Nothing
   else
@@ -70,10 +83,10 @@ function rand(t::Union{Type{DAMap},Type{TPSAMap}}; spin::Union{Bool,Nothing}=not
     end
   end
 
-  if isnothing(stochastic)
+  if isnothing(FD)
     V = Nothing
   else
-    if stochastic
+    if FD
       V = Matrix{Float64}
     else
       V = Nothing
@@ -161,15 +174,15 @@ function read_fpp_map(file; idpt::Union{Nothing,Bool}=nothing)
 
   # Check if map has stochasticity
   stoch_idx = findfirst(t->t=="Stochastic", data)
-  if !isnothing(stoch_idx) && stoch_idx[2]== 1 # stochastic in first column meaning no "No"
-    stochastic = true
+  if !isnothing(stoch_idx) && stoch_idx[2]== 1 # FD in first column meaning no "No"
+    FD = true
   else
-    stochastic=nothing
+    FD=nothing
   end
 
   # Make the TPSA
   d = Descriptor(nv, no, np, no)
-  m = complex(DAMap(use=d,idpt=idpt,stochastic=stochastic)) #repeat([ComplexTPS(use=d)], nv), Q=Quaternion([ComplexTPS(1,use=d), repeat([ComplexTPS(use=d)], 3)...]))
+  m = complex(DAMap(use=d,idpt=idpt,FD=FD)) #repeat([ComplexTPS(use=d)], nv), Q=Quaternion([ComplexTPS(1,use=d), repeat([ComplexTPS(use=d)], 3)...]))
 
   idx=3
   data=data[3:end,:]
@@ -198,8 +211,8 @@ function read_fpp_map(file; idpt::Union{Nothing,Bool}=nothing)
   # dont forget params
   m.x[nv+1:nn] .= complexparams(d)
 
-  # stochastic
-  if !isnothing(stochastic)
+  # FD
+  if !isnothing(FD)
     idx = findfirst(t->t=="Stochastic", data)[1]
     data = data[idx+1:end,:]
     for i=1:size(data, 1)
diff --git a/src/utils/quaternion.jl b/src/utils/quaternion.jl
index 3633582..e8784e7 100644
--- a/src/utils/quaternion.jl
+++ b/src/utils/quaternion.jl
@@ -13,6 +13,8 @@ Quaternion(Q::Quaternion{T}) where T <: Number = Quaternion(map(x->T(x), Q.q))
 Quaternion(t::T) where T <: Number = Quaternion([one(t), zero(t), zero(t), zero(t)])
 Quaternion(::Nothing) = Quaternion{Nothing}(Nothing[])
 
+eltype(::Type{Quaternion{T}}) where T = T
+
 ==(Q1::Quaternion, Q2::Quaternion) = Q1.q == Q2.q
 
 """
diff --git a/src/vectorfield/ctors.jl b/src/vectorfield/ctors.jl
index 3fde5fc..4197945 100644
--- a/src/vectorfield/ctors.jl
+++ b/src/vectorfield/ctors.jl
@@ -66,11 +66,11 @@ function VectorField{T,U}(h::T; Q::Union{U,Nothing}=nothing, work_low::Vector{<:
 end
 
 """
-    VectorField{T,U}(u::UndefInitializer; use::UseType=GTPSA.desc_current) where {T,U}
+    VectorField{T,U}(u::UndefInitializer; use=GTPSA.desc_current) where {T,U}
 
 Creates an undefined `VectorField{T,U}` with same number of variables as `use`.
 """
-function VectorField{T,U}(u::UndefInitializer; use::UseType=GTPSA.desc_current) where {T,U}
+function VectorField{T,U}(u::UndefInitializer; use=GTPSA.desc_current) where {T,U}
   desc = getdesc(use)
   nv = numvars(desc)
   x = Vector{T}(undef, nv)
@@ -137,7 +137,7 @@ function zero(F::VectorField{T,U}) where {T,U}
   return VectorField(x,Q)
 end
 
-function zero(::Type{VectorField{T,U}}; use::UseType=GTPSA.desc_current) where {T,U}
+function zero(::Type{VectorField{T,U}}; use=GTPSA.desc_current) where {T,U}
   desc = getdesc(use)
   nn = numnn(desc)
   nv = numvars(desc)
diff --git a/src/vectorfield/methods.jl b/src/vectorfield/methods.jl
index 629b2f3..63db286 100644
--- a/src/vectorfield/methods.jl
+++ b/src/vectorfield/methods.jl
@@ -1,5 +1,7 @@
 # --- copy! ---
-function copy!(F::VectorField{T,U}, F1::Union{VectorField{T,U},TaylorMap{S,T,U,V}}) where {S,T,U,V}
+function copy!(F::VectorField, F1::Union{VectorField,TaylorMap})
+  checkinplace(F, F1)
+
   desc = getdesc(F)
   nv = numvars(desc)
 
@@ -17,7 +19,7 @@ function copy!(F::VectorField{T,U}, F1::Union{VectorField{T,U},TaylorMap{S,T,U,V
 end
 
 # --- complex ---
-function complex(F::VectorField{T,U}) where {T,U}
+function complex(F::VectorField)
   desc = getdesc(F)
   nn = numnn(desc)
   nv = numvars(desc)
@@ -51,7 +53,7 @@ liebra!(na::Cint, m1::Vector{Ptr{RTPSA}}, m2::Vector{Ptr{RTPSA}}, m3::Vector{Ptr
 liebra!(na::Cint, m1::Vector{Ptr{CTPSA}}, m2::Vector{Ptr{CTPSA}}, m3::Vector{Ptr{CTPSA}}) = (@inline; GTPSA.mad_ctpsa_liebra!(na, m1, m2, m3))
 
 """
-    lb(F::VectorField{T,U}, H::VectorField{T,U}) where {T,U}
+    lb(F::VectorField, H::VectorField)
 
 Computes the Lie bracket of the vector fields `F` and `H`. Explicitly, and including 
 spin (with the lower case letter for the quaternion of the vector field), this is:
@@ -60,7 +62,9 @@ spin (with the lower case letter for the quaternion of the vector field), this i
 
 where `[h,f] = h*f - f*h` is just a quaternion commutator. See Equation 44.52 in the Bmad manual
 """
-function lb(F::VectorField{T,U}, H::VectorField{T,U}) where {T,U}
+function lb(F::VectorField, H::VectorField)
+  checkop(F,H)
+
   G = zero(F)
   lb!(G,F,H)
   return G
@@ -68,7 +72,7 @@ end
 
 
 """
-    lb!(G::VectorField{T,U}, F::VectorField{T,U}, H::VectorField{T,U}; work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_vf_work_low(F), work_Q::Union{Nothing,U}=nothing) where {T,U}
+    lb!(G::VectorField, F::VectorField, H::VectorField; work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_vf_work_low(F), work_Q::Union{Quaternion,Nothing}=nothing)
 
 Computes the Lie bracket of the vector fields `F` and `H`, and stores the result in `G`. 
 See `lb` for more details.
@@ -77,7 +81,11 @@ See `lb` for more details.
 - `work_low` -- Tuple of 3 vectors of type `lowtype(T)` of length `>=nv`
 - `work_Q`   -- `Quaternion{T}` if spin is included in the vector field, else `nothing`
 """
-function lb!(G::VectorField{T,U}, F::VectorField{T,U}, H::VectorField{T,U}; work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_lb_work_low(F), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {T,U}
+function lb!(G::VectorField, F::VectorField, H::VectorField; work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_lb_work_low(F), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F))
+  checkinplace(G, F, H)
+
+  T = eltype(G.x)
+
   nv = numvars(F)
   Fx_low = work_low[1]
   Hx_low = work_low[2]
@@ -125,22 +133,23 @@ end
 # --- Lie operator (VectorField) acting on DAMap ---
 # GTPSA provides fgrad with is used, but nothing for spin
 """
-    *(F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}) where {S,T,U,V,W}
+    *(F::VectorField, m1::DAMap)
 
 Calculates a Lie operator `VectorField` acting on a `DAMap`. Explicity, if spin is 
 included that is `F * m = (F.x, F.Q) * (m.x, m.Q) = (F.x ⋅ ∇ m.x , F.x ⋅ ∇ m.Q + m.Q*F.Q)`
 """
-function *(F::VectorField{T,U}, m1::Union{DAMap{S,T,U,V,W},UniformScaling}) where {S,T,U,V,W}
+function *(F::VectorField, m1::Union{DAMap,UniformScaling})
   if m1 isa UniformScaling 
-    m1 = one(DAMap{numtype(T),T,U,Nothing,Nothing},use=getdesc(F))
+    m1 = one(DAMap{numtype(eltype(F.x)),eltype(F.x),typeof(F.Q),Nothing,Nothing},use=getdesc(F))
   end
+  checkop(F, m1)
   m = zero(m1)
   mul!(m, F, m1)
   return m
 end
 
 """
-    mul!(m::DAMap{S,T,U,V,W}, F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}; work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {S,T,U,V,W}
+    mul!(m::DAMap, F::VectorField, m1::DAMap; work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F))
 
 Computes the Lie operator `F` acting on a `DAMap` `m1`, and stores the result in `m`.
 Explicity, that is `F * m = (F.x, F.Q) * (m.x, m.Q) = (F.x ⋅ ∇ m.x , F.x ⋅ ∇ m.Q + m.Q*F.Q)`
@@ -149,9 +158,10 @@ Explicity, that is `F * m = (F.x, F.Q) * (m.x, m.Q) = (F.x ⋅ ∇ m.x , F.x ⋅
 - `work_low` -- Vector with eltype `lowtype(T)` and length `>=nv`
 - `work_Q`   -- `Quaternion{T}` if spin is included in the vector field, else `nothing`
 """
-function mul!(m::DAMap{S,T,U,V,W}, F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}; work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {S,T,U,V,W}
-  checkidpt(m,m1)
-  
+function mul!(m::DAMap, F::VectorField, m1::DAMap; work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F))
+  checkinplace(m, F, m1)
+
+  T = eltype(m.x)
   nv = numvars(F)
 
   @assert length(work_low) >= nv "Incorrect length for work_low; received $(length(work_low)), should be >=$nv"
@@ -188,10 +198,12 @@ end
 # --- exp(F)*m ---
 # While GTPSA provides this for only the orbital part, it separately converges 
 # each variable and does not include spin. Here I include spin and do the entire map at once.
-function exp(F::VectorField{T,U}, m1::Union{UniformScaling,DAMap{S,T,U,V,W}}=I) where {S,T,U,V,W}
+function exp(F::VectorField, m1::Union{UniformScaling,DAMap}=I)
   if m1 isa UniformScaling 
-    m1 = one(DAMap{numtype(T),T,U,Nothing,Nothing},use=getdesc(F))
+    m1 = one(DAMap{numtype(eltype(F.x)),eltype(F.x),typeof(F.Q),Nothing,Nothing},use=getdesc(F))
   end
+  checkop(F,m1)
+
   m = zero(m1)
   exp!(m, F, m1)
   return m
@@ -199,7 +211,7 @@ end
 
 
 """
-    exp!(m::DAMap{S,T,U,V,W}, F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}; work_maps::Tuple{Vararg{DAMap{S,T,U,V,W}}}=(zero(m1),zero(m1)), work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {S,T,U,V,W}
+    exp!(m::DAMap, F::VectorField, m1::DAMap; work_maps::Tuple{Vararg{DAMap}}=(zero(m1),zero(m1)), work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F))
 
 Computes `exp(F)*m1`, and stores the result in `m`. Explicity, this is
 `exp(F)*m1 = m1 + F*m1 + 1/2*F*(F*m1) + 1/6*F*(F*(F*m1)) + ...`, where `*` is
@@ -212,9 +224,10 @@ number of iterations is equal to the number of higher orders left.
 - `work_low` -- Vector with eltype `lowtype(T)` and length `>=nv`
 - `work_Q`   -- `Quaternion{T}` if spin is included in the vector field, else `nothing`
 """
-function exp!(m::DAMap{S,T,U,V,W}, F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}; work_maps::Tuple{Vararg{DAMap{S,T,U,V,W}}}=(zero(m1),zero(m1)), work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {S,T,U,V,W}
-  checkidpt(m, m1, work_maps...)
-  
+function exp!(m::DAMap, F::VectorField, m1::DAMap; work_maps::Tuple{Vararg{DAMap}}=(zero(m1),zero(m1)), work_low::Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}=Vector{lowtype(first(F.x))}(undef, numvars(F)), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F))
+  checkinplace(m, F, m1, work_maps...)
+
+  T = eltype(m.x)
   nv = numvars(F)
   #GTPSA.mad_ctpsa_exppb!(nv, map(t->t.tpsa, F.x), map(t->t.tpsa, view(m1.x, 1:nv)), map(t->t.tpsa, view(m.x,1:nv)))
   #return
@@ -272,7 +285,7 @@ end
 # 2 work_vfs
 
 """
-    log!(F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}; work::Tuple{DAMap{S,T,U,V,W},DAMap{S,T,U,V,W},DAMap{S,T,U,V,W},VectorField{T,U},VectorField{T,U}}=prep_log_work(m1), work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_lb_work_low(F), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {S,T,U,V,W}
+    log!(F::VectorField, m1::DAMap; work::Tuple{DAMap,DAMap,DAMap,VectorField,VectorField}=prep_log_work(m1), work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_lb_work_low(F), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F)) 
 
 Computes the log of the map `m1` - that is, calculates the `VectorField` `F` that 
 would represent the map `m1` as a Lie exponent `exp(F)` - and stores the result in `F`.
@@ -283,8 +296,8 @@ The map `m1` should be close to the identity for this to converge quickly.
 - `work_low` -- Vector with eltype `lowtype(T)` and length `>=nv`
 - `work_Q`   -- `Quaternion{T}` if spin is included in the vector field, else `nothing`
 """
-function log!(F::VectorField{T,U}, m1::DAMap{S,T,U,V,W}; work::Tuple{DAMap{S,T,U,V,W},DAMap{S,T,U,V,W},DAMap{S,T,U,V,W},VectorField{T,U},VectorField{T,U}}=prep_log_work(m1), work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_lb_work_low(F), work_Q::Union{U,Nothing}=prep_vf_work_Q(F)) where {S,T,U,V,W}
-  checkidpt(m1, work[1:3]...)
+function log!(F::VectorField, m1::DAMap; work::Tuple{DAMap,DAMap,DAMap,VectorField,VectorField}=prep_log_work(m1), work_low::Tuple{Vararg{Vector{<:Union{Ptr{RTPSA},Ptr{CTPSA}}}}}=prep_lb_work_low(F), work_Q::Union{Quaternion,Nothing}=prep_vf_work_Q(F))
+  checkinplace(F, m1, work...)
   
   nv = numvars(m1)
   nmax = 100
@@ -368,14 +381,14 @@ end
 
 
 """
-    log(m1::DAMap{S,T,U,V}) where {S,T,U,V}
+    log(m1::DAMap)
 
 Computes the log of the map `m1` - that is, calculates the `VectorField` `F` that 
 would represent the map `m1` as a Lie exponent `exp(F)`. The map `m1` should be close 
 to the identity for this to converge quickly.
 """
-function log(m1::DAMap{S,T,U,V}) where {S,T,U,V}
-  F = zero(VectorField{T,U},use=m1)
+function log(m1::DAMap)
+  F = zero(VectorField{eltype(m1.x),typeof(m1.Q)},use=m1)
   log!(F,m1)
   return F
 end
diff --git a/src/vectorfield/operators.jl b/src/vectorfield/operators.jl
index 355fefe..9938f19 100644
--- a/src/vectorfield/operators.jl
+++ b/src/vectorfield/operators.jl
@@ -16,7 +16,7 @@ for ops = (("add!", :+), ("sub!",:-))
 @eval begin
   
 function $(Meta.parse(ops[1]))(F::VectorField, F1::Union{VectorField,TaylorMap}, F2::Union{VectorField,TaylorMap}; dospin::Bool=true)
-  checkspin(F,F1,F2)
+  checkinplace(F, F1, F2)
   
   nv = numvars(F)
 
@@ -34,7 +34,7 @@ function $(Meta.parse(ops[1]))(F::VectorField, F1::Union{VectorField,TaylorMap},
 end
 
 function $(Meta.parse(ops[1]))(F::VectorField, J::UniformScaling, F1::Union{VectorField,TaylorMap}; dospin::Bool=true)
-  checkspin(F,F1)
+  checkinplace(F, F1)
   
   nv = numvars(F)
 
@@ -54,7 +54,7 @@ function $(Meta.parse(ops[1]))(F::VectorField, J::UniformScaling, F1::Union{Vect
 end
 
 function $(Meta.parse(ops[1]))(F::VectorField, F1::Union{VectorField,TaylorMap}, J::UniformScaling; dospin::Bool=true)
-  checkspin(F,F1)
+  checkinplace(F, F1)
   
   nv = numvars(F)
 
@@ -74,7 +74,8 @@ function $(Meta.parse(ops[1]))(F::VectorField, F1::Union{VectorField,TaylorMap},
 end
 
 function $(ops[2])(F1::VectorField, F2::VectorField)
-  checkspin(F1,F1)
+  checkop(F1, F2)
+
   # Promote if necessary:
   if eltype(F1.x) == ComplexTPS
     F = zero(F1)
@@ -107,7 +108,7 @@ end
 for ops = (("add!", :+), ("sub!",:-), ("mul!",:*), ("div!",:/))
 @eval begin
 function $(Meta.parse(ops[1]))(F::VectorField, a::Number, F1::VectorField; dospin::Bool=true)
-  checkspin(F,F1)
+  checkinplace(F, F1, a)
   
   nv = numvars(F)
 
@@ -124,7 +125,8 @@ function $(Meta.parse(ops[1]))(F::VectorField, a::Number, F1::VectorField; dospi
 end
 
 function $(Meta.parse(ops[1]))(F::VectorField, F1::VectorField, a::Number; dospin::Bool=true)
-  checkspin(F,F1)
+  checkinplace(F, F1, a)
+
   nv = numvars(F)
 
   for i=1:nv
diff --git a/test/runtests.jl b/test/runtests.jl
index c9297bc..6b61a83 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,20 +4,20 @@ using Test
 @testset "NonlinearNormalForm.jl" begin
     d = Descriptor(1,2)
     x1 = vars()[1]
-    m1 = DAMap(x=[1+2*x1+2*x1^2], x0=[4])
-    m2 = DAMap(x=[1+2*x1+2*x1^2], x0=[3])
+    m1 = DAMap(x=[1+2*x1+2*x1^2], x0=[4.])
+    m2 = DAMap(x=[1+2*x1+2*x1^2], x0=[3.])
 
     mt1 = TPSAMap(m1)
     mt2 = TPSAMap(m2)
     
     tol = 1e-10
 
-    @test norm(m2∘m1 - DAMap(x=[1+4*x1+12*x1^2], x0=[4])) < tol
-    @test norm(mt2∘mt1 - TPSAMap(x=[5-12*x1-4*x1^2], x0=[4])) < tol
+    @test norm(m2∘m1 - DAMap(x=[1+4*x1+12*x1^2], x0=[4.])) < tol
+    @test norm(mt2∘mt1 - TPSAMap(x=[5-12*x1-4*x1^2], x0=[4.])) < tol
     @test norm(m2^3 - m2∘m2∘m2) < tol
     @test norm(mt2^3 - mt2∘mt2∘mt2) < tol
-    @test norm(m2^3 - DAMap(x=[1+8*x1+56*x1^2], x0=[3])) < tol
-    @test norm(mt2^3 - TPSAMap(x=[13+8*x1-8*x1^2], x0=[3])) < tol
+    @test norm(m2^3 - DAMap(x=[1+8*x1+56*x1^2], x0=[3.])) < tol
+    @test norm(mt2^3 - TPSAMap(x=[13+8*x1-8*x1^2], x0=[3.])) < tol
 
     # with temporary inverter"
     #@test norm(m1^-3 - DAMap(x=[4+0.125*x1-0.109375*x1^2],x0=[1])) < tol