text:
export MGsetup,transposeHierarchy,adjustMemoryForNumRHS,getSPAIprec,getFWInterp,replaceMatrixInHierarchy
export restrictCellCenteredVariables,restrictNodalVariables2;
single = false;
function MGsetup(ATf::Union{SparseMatrixCSC,multilevelOperatorConstructor},Mesh::RegularMesh,param::MGparam,rhsType::DataType = Float64,nrhs::Int64 = 1,verbose::Bool=false)
Ps      	= Array{SparseMatrixCSC}(param.levels-1);
Rs      	= Array{SparseMatrixCSC}(param.levels-1);
As 			= Array{SparseMatrixCSC}(param.levels);
Meshes  	= Array{RegularMesh}(param.levels);
relaxPrecs 	= Array{Any}(param.levels);

PDEparam = 0;
if isa(ATf,SparseMatrixCSC)==true
	As[1] = ATf;
else
	As[1] = ATf.getOperator(Mesh,ATf.param)';
	PDEparam = ATf.param;
end
Meshes[1] = Mesh;
n = Mesh.n;
AT = As[1];
Cop = nnz(AT);
for l = 1:(param.levels-1)
	if verbose
		tic()
	end
    AT = As[l];

	PT = 0;
	RT = 0;
	withCellsBlock = false;
	if param.transferOperatorType=="SystemsFacesMixedLinear"
		withCellsBlock = true;
	end

	if param.transferOperatorType=="FullWeighting"
		geometric = !isa(ATf,SparseMatrixCSC);
		(P,nc) = getFWInterp(n+1,geometric);
		if isa(ATf,SparseMatrixCSC)
			RT = P;
			PT = P';
		else
			RT = (0.5^Meshes[l].dim)*P;
			PT = P';
		end
		nc = nc-1;
		P = 0;
	elseif param.transferOperatorType=="SystemsFacesLinear" || param.transferOperatorType=="SystemsFacesMixedLinear"
		(P,R,nc) = getLinearOperatorsSystemsFaces(n,withCellsBlock);
		PT = P';
		RT = R';
		P = 0;
		R = 0;
	end

	if param.relaxType=="Jac" || param.relaxType=="Jac-GMRES"
		d = param.relaxParam./diag(AT);
		d = spdiagm(d);
		if single
			relaxPrecs[l] = convert(SparseMatrixCSC{Complex64,Int32},d);# here we need to take the conjugate for the SpMatVec, but we give At instead of A so it cancels
		else
			relaxPrecs[l] = d;
		end

	elseif param.relaxType=="SPAI"
		relaxPrecs[l] = spdiagm(param.relaxParam*getSPAIprec(AT)); # here we need to take the conjugate for the SpMatVec, but we give At instead of A so it cancels
	elseif param.relaxType=="VankaFaces"
		relaxPrecs[l] = getVankaFacesPreconditioner(AT,Meshes[l],param.relaxParam,withCellsBlock);
	else
		error("Unknown relaxation type !!!!");
	end

	if (size(PT,1)==size(PT,2))
		if verbose; println(string("Stopped Coarsening at level ",l)); end
		param.levels = l;
		As = As[1:l];
		Ps = Ps[1:l-1];
		Rs = Rs[1:l-1];
		Meshes = Meshes[1:l];
		relaxPrecs = relaxPrecs[1:l];
		break;
	else
		Ps[l] = PT;
		Rs[l] = RT;
		Meshes[l+1] = getRegularMesh(Meshes[l].domain,nc);
		if isa(ATf,SparseMatrixCSC)
			# if param.transferOperatorType=="SystemsFacesLinear"
				# Rinj = getInjectionOperatorsSystemsFaces(n,withCellsBlock)';
				# Act = Ps[l]*AT*Rinj;
			# else
				Act = Ps[l]*AT*Rs[l];
			# end
		else
			PDEparam = ATf.restrictParams(Meshes[l],Meshes[l+1],PDEparam,l);
			Act = ATf.getOperator(Meshes[l+1],PDEparam)';
		end
		As[l+1] = Act;
		Cop = Cop + nnz(Act);
		if verbose; println("MG setup: ",n," cells took:",toq()); end;
		n = nc;
	end
end
if verbose
	tic()
	println("MG setup: Operator complexity = ",Cop/nnz(As[1]));
end

param.As = As;
param.Meshes = Meshes;

defineCoarsestAinv(param,As[end]);

if verbose
	println("MG setup coarsest ",param.coarseSolveType,": ",n," cells , done LU in ",toq());
end
param.Ps = Ps;
param.Rs = Rs;
param.relaxPrecs = relaxPrecs;
param = adjustMemoryForNumRHS(param,rhsType,nrhs,verbose);
param.doTranspose = 0;
return;
end

function adjustMemoryForNumRHS(param::MGparam,rhsType::DataType = Float64,nrhs::Int64 = 1,verbose::Bool=false)
if length(param.As)==0
	error("The Hierarchy is empty - run a setup first.")
end
good = true;
if length(param.memCycle) > 0
	if size(param.memCycle[1].x,2) != nrhs
		good = false
	elseif param.relaxType == "Jac-GMRES"
		if length(param.memRelax)==0
			good = false;
		end
	elseif param.cycleType == 'K'
		if length(param.memRelax)==0
			good = false;
		end
	end
else
	good = false;
end
if good
	return param;
end

if param.cycleType == 'K'
	memKcycle = Array{FGMRESmem}(max(param.levels-2,0));
else
	memKcycle = Array{FGMRESmem}(0);
end

if param.relaxType == "Jac-GMRES"
	memRelax = Array{FGMRESmem}(param.levels-1);
else
	memRelax = Array{FGMRESmem}(0);
end

memCycle = Array{CYCLEmem}(param.levels);

param.memRelax  = memRelax;
param.memKcycle = memKcycle;
param.memCycle  = memCycle;

N = size(param.As[1],2);

for l = 1:(param.levels-1)
	N = size(param.As[l],2);
	needZ = param.relaxType=="Jac-GMRES";
	memCycle[l] = getCYCLEmem(N,nrhs,rhsType,true);
	if param.relaxType=="Jac-GMRES"
		maxRelax = max(param.relaxPre(l),param.relaxPost(l));
		if nrhs == 1
			memRelax[l] = getFGMRESmem(N,true,rhsType,maxRelax,1);
		else
			# memRelax[l] = getBlockFGMRESmem(N,nrhs,false,rhsType,maxRelax);
			# memRelax[l] = getFGMRESmem(N,false,rhsType,maxRelax,nrhs);
			memRelax[l] = getFGMRESmem(N,true,rhsType,maxRelax,1);
		end
	end

	if l > 1 && param.cycleType=='K'
		if nrhs == 1
			memKcycle[l-1] = getFGMRESmem(N,true,rhsType,2,1);
		else
			# memKcycle[l-1] = getBlockFGMRESmem(N,nrhs,true,rhsType,2);
			memKcycle[l-1] = getFGMRESmem(N,true,rhsType,2,nrhs);
		end

	end
end
memCycle[end] = getCYCLEmem(size(param.As[end],2),nrhs,rhsType,false); # no need for residual on the coarsest level...
memCycle[end].b = memCycle[end].r;
param.memRelax = memRelax;
param.memKcycle = memKcycle;
param.memCycle = memCycle;
return param;
end

function replaceMatrixInHierarchy(param,AT::SparseMatrixCSC,verbose::Bool=false)

param.As[1] = AT;
Cop = nnz(AT);
for l = 1:(param.levels-1)
	if verbose
		tic()
	end
    AT = param.As[l];

	if param.relaxType=="Jac" || param.relaxType=="Jac-GMRES"
		d = param.relaxParam./diag(AT);
		d = spdiagm(d);
		param.relaxPrecs[l] = d;
		if single
			param.relaxPrecs[l] = convert(SparseMatrixCSC{Complex64,Int32},d);# here we need to take the conjugate for the SpMatVec, but we give At instead of A so it cancels
		end

	elseif param.relaxType=="SPAI"
		d = spdiagm(d);
	elseif param.relaxType=="VankaFaces"
		param.relaxPrecs[l] = getVankaFacesPreconditioner(AT,Meshes[l],param.relaxParam,withCellsBlock);
	else
		error("Unknown relaxation type !!!!");
	end
	Act = param.Ps[l]*AT*param.Rs[l];
	param.As[l+1] = Act;
	Cop = Cop + nnz(Act);
	if verbose; println("MG setup: took:",toq()); end;
end
if verbose
	tic()
	println("MG setup: Operator complexity = ",Cop/nnz(param.As[1]));
end

defineCoarsestAinv(param,param.As[end]);

if verbose
	println("MG setup coarsest ",param.coarseSolveType,":  done LU in ",toq());
end
param.doTranspose = 0;
return;
end

function transposeHierarchy(param::MGparam,verbose::Bool=false)
# This is a shortened MG setup. It assumes that a setup was previously applied
# to the matrix in param.As[1] (which is held transposed). It keeps the same Ps and memory allocations.
# Now all hierarchy is transposed. doTranspose is set to 1. MUMPS do not replace the factorization cause it knows to solve transpose.

param.As[1] = param.As[1]';
param.doTranspose = mod(param.doTranspose+1,2);

for l = 1:(param.levels-1)
	if verbose
		println("MG transpose level: ",l);
		tic()
	end

    if param.relaxType=="Jac" || param.relaxType=="Jac-GMRES" || param.relaxType=="SPAI"
		param.relaxPrecs[l] = conj(param.relaxPrecs[l]);
	else
		error("Not supported");
	end
	R = param.Ps[l]';
	param.Ps[l] = param.Rs[l]';
	param.Rs[l] = R;

	param.As[l+1]  = param.As[l+1]';
end
if verbose
	tic()
end
if param.coarseSolveType == "MUMPS"
	# when applying mumps, no need to transpose the factorization
elseif param.coarseSolveType == "BiCGSTAB" || param.coarseSolveType == "GMRES"
	param.LU = conj(param.LU);
else
	destroyCoarsestLU(param);
	param.LU = lufact(param.As[end]');
end
if verbose
	println("MG transpose hierarchy: done LU");
	toc();
end
return;
end

export defineCoarsestAinv;
function defineCoarsestAinv(param::MGparam,AT::SparseMatrixCSC)
if param.coarseSolveType == "MUMPS"
	param.LU = factorMUMPS(AT',0,0);
elseif param.coarseSolveType == "BiCGSTAB" || param.coarseSolveType == "GMRES"
	d = spdiagm(param.relaxParam./diag(AT));
	param.LU = d;
	if single
		param.LU = convert(SparseMatrixCSC{Complex64,Int32},d);# here we need to take the conjugate for the SpMatVec, but we give At instead of A so it cancels
	end

# elseif param.coarseSolveType[1:7] == "DDNodal"
	# DDNumCells = parseDD(param.coarseSolveType,param.Meshes[end].dim);
	# param.LU,nnzDD = setupDD(AT,param.Meshes[end],getNodalIndicesOfCell,DDNumCells,1e-5);
# elseif param.coarseSolveType[1:7] == "DDFaces"
	# DDNumCells = parseDD(param.coarseSolveType,param.Meshes[end].dim);
	# param.LU,nnzDD = setupDD(AT,param.Meshes[end],getFacesStaggeredIndicesOfCell, DDNumCells,1e-5);
	# # println("Cop coarsest (DomainDecomp): ",nnzDD./nnz(param.As[1]));
else
	param.LU = lufact(AT');
	# nnzLU = nnz(param.LU[:L]) + nnz(param.LU[:U]);
	# println("Cop coarsest (Julia): ",nnzLU./nnz(param.As[1]));
end
end

# function parseDD(DDstring::String,dim::Int64)
# if dim==2
	# DDNumCells = [parse(Int64,DDstring[end-1]),parse(Int64,DDstring[end])];
# else
	# DDNumCells = [parse(Int64,DDstring[end-2]),parse(Int64,DDstring[end-1]),parse(Int64,DDstring[end])];
# end
# return DDNumCells;
# end

# Bilinear "Full Weighting" prolongation
function getFWInterp(n_nodes::Array{Int64,1},geometric::Bool=false)
# n here is the number of NODES!!!
(P1,nc1) = get1DFWInterp(n_nodes[1],geometric);
(P2,nc2) = get1DFWInterp(n_nodes[2],geometric);
if length(n_nodes)==3
	(P3,nc3) = get1DFWInterp(n_nodes[3],geometric);
end
if length(n_nodes)==2
	P = kron(P2,P1);
	nc = [nc1,nc2];
else
	P = kron(P3,kron(P2,P1));
	nc = [nc1,nc2,nc3];
end
return P,nc
end

function get1DFWInterp(n_nodes::Int64,geometric)
# n here is the number of NODES!!!
oddDim = mod(n_nodes,2);
if n_nodes > 2
	halfVec = 0.5*ones(n_nodes-1);
	P = spdiagm((halfVec,ones(n_nodes),halfVec),[-1,0,1],n_nodes,n_nodes);
    if oddDim == 1
        P = P[:,1:2:end];
    else
		if geometric
			P = speye(n_nodes);
			warn("getFWInterp: in geometric mode we stop coarsening because num cells does not divide by two");
		else
			P = P[:,[1:2:end;end]];
			P[end-1:end,end-1:end] = speye(2);
#         	P = P[:,1:2:end];
#         	P[end,end-1:end] = [-0.5,1.5];
		end
    end
else
    P = speye(n_nodes);
end
nc = size(P,2);
if single
	P = convert(SparseMatrixCSC{Float32,Int32},P);
end
return P,nc
end

function getSPAIprec(A::SparseMatrixCSC)
s = vec(sum(real(A).^2,2) + sum(imag(A).^2,2));
Q = conj(diag(A))./s;
end

########################### GEOMETRIC MULTIGRID STUFF ###############################

function restrictCellCenteredVariables(rho::ArrayTypes,n::Array{Int64})
R,nc = getRestrictionCellCentered(n);
rho_c = R*rho[:];
return rho_c;
## TODO: make this more efficient...
end

export restrictNodalVariables
function restrictNodalVariables(rho::ArrayTypes,n_nodes::Array{Int64})
P,nc = getFWInterp(n_nodes,true);
rho_c = zeros(eltype(rho),size(P,2));
rho_c[:] = (0.5^length(n_nodes))*(P'*rho[:]);
return rho_c;
## TODO: make this more efficient...
end

function restrictNodalVariables2(rho::ArrayTypes,n_nodes::Array{Int64})
P,nc = getFWInterp(n_nodes,true);
R1,nc1 = get1DNodeFullWeightRestriction(n_nodes[1]-1);
R2,nc2 = get1DNodeFullWeightRestriction(n_nodes[2]-1);
R = kron(R2,R1);
# rho_c = zeros(eltype(rho),size(P,2));
# println(size(R))
# println(size(rho[:]))

rho_c = R*rho[:];
return rho_c;
## TODO: make this more efficient...
end


text:
# Basic UserGroupInformation to store user information and tokens,
# realUser: what it is at client side
# effectiveUser: set this to specify what this should be proxied as at the server end
# ref: https://www.opencore.com/blog/2016/5/user-name-handling-in-hadoop/
mutable struct UserGroupInformation
    userinfo::UserInformationProto
    credentials::Credentials

    function UserGroupInformation(username::AbstractString=default_username(); proxy::Bool=false, proxyuser::AbstractString=username)
        userinfo = proxy ? UserInformationProto(realUser=username, effectiveUser=proxyuser) : UserInformationProto(realUser=username)
        ugi = new(userinfo, Credentials())
        haskey(ENV, "HADOOP_TOKEN_FILE_LOCATION") && read_credentials!(ENV["HADOOP_TOKEN_FILE_LOCATION"]; credentials=ugi.credentials)
        ugi
    end
end

function default_username()
    ("HADOOP_USER_NAME" in keys(ENV)) && (return ENV["HADOOP_USER_NAME"])
    ("USER" in keys(ENV)) && (return ENV["USER"])
    error("Can not determine user information. Either HADOOP_USER_NAME or USER must be set.")
end

add_token!(ugi::UserGroupInformation, token::TokenProto) = add_token!(ugi.credentials, token)
add_token!(ugi::UserGroupInformation, alias::AbstractString, token::TokenProto) = add_token!(ugi.credentials, alias, token)

username(userinfo::UserInformationProto) = hasproperty(userinfo, :realUser) ? userinfo.realUser : userinfo.effectiveUser
username(ugi::UserGroupInformation) = username(ugi.userinfo)

find_tokens(ugi::UserGroupInformation; alias::AbstractString="", kind::AbstractString="") = find_tokens(ugi.credentials; alias=alias, kind=kind)

function show(io::IO, ugi::UserGroupInformation)
    uinfo = ugi.userinfo
    print(io, "User:")
    hasproperty(uinfo, :realUser) && print(io, ' ', uinfo.realUser)
    hasproperty(uinfo, :effectiveUser) && print(io, " (", uinfo.effectiveUser, ')')
    isempty(ugi.credentials.tokens) || print(io, " with ", length(ugi.credentials.tokens), " tokens")
    nothing
end

ext: jl


text:
# * Field arithmetic [1/1]
#   For calculations, a time-base is necessary ("normalized time"), with
#   respect to which all harmonic motions &c are analyzed.  When
#   combining fields of commensurate frequencies, it is easy to
#   establish a common time-base, most likely the fundamental wavelength
#   will be the obvious choice -- such as when adding an IR field and its
#   harmonic components as generated through e.g. HHG:
#
#   \[E(t) = \sum_q E_{2q + 1}(t)\sin[(2q+1)\omega t].\]
#
#   However, when adding incommensurate frequencies, there is no obvious
#   choice, so the user has to specify the time-base explicitly -- maybe
#   by the order in which the fields are added?
#
#   There should also be some helper routines that, when discretizing
#   the time axis, estimate whether all harmonic components of interest
#   are satisfactorily resolved.
#
# ** Sum fields

import Base: +

"""
    SumField(a, b)

The linear combination of two fields `a` and `b`.

# Example

```jldoctest
julia> @field(A) do
           I₀ = 1.0
           T = 2.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> @field(B) do
           I₀ = 0.5
           T = 1.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 5.0000e-01 au = 1.7547226e16 W cm⁻² =>
    - E₀ = 7.0711e-01 au = 363.6089 GV m⁻¹
    - A₀ = 0.1125 au
  – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±1.00σ)
  – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm

julia> A+B
┌ Linearly polarized field with
│   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
│     - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
│     - A₀ = 0.3183 au
│   – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
│   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
│   – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
⊕
│ Linearly polarized field with
│   - I₀ = 5.0000e-01 au = 1.7547226e16 W cm⁻² =>
│     - E₀ = 7.0711e-01 au = 363.6089 GV m⁻¹
│     - A₀ = 0.1125 au
│   – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV)
│   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±1.00σ)
└   – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm
```
"""
struct SumField{A,B} <: AbstractField
    a::A
    b::B
    function SumField(a::A, b::B) where {A<:AbstractField, B<:AbstractField}
        da = dimensions(a)
        db = dimensions(b)
        da == db ||
            throw(ArgumentError("Cannot add fields of different dimensionality, got $(da) and $(db)"))
        new{A,B}(a, b)
    end
end

function show(io::IO, f::SumField)
    a_str = split(string(f.a), "\n")
    b_str = split(string(f.b), "\n")

    for (s,l) in zip("┌" * repeat("│", length(a_str)-1), a_str)
        write(io, "$s $l\n")
    end

    write(io, "⊕\n")

    for (s,l) in zip(repeat("│", length(b_str)-1) * "└", b_str)
        write(io, "$s $l\n")
    end
end

function +(a::AbstractField, b::AbstractField)
    polarization(a) == polarization(b) ||
        throw(ArgumentError("Cannot add fields of different polarization"))
    SumField(a, b)
end

vector_potential(f::SumField, t::Number) =
    vector_potential(f.a, t) + vector_potential(f.b, t)

polarization(f::SumField) = polarization(f.a)

function span(f::SumField)
    sa = span(f.a)
    sb = span(f.b)
    min(sa.left,sb.left)..max(sa.right,sb.right)
end
steps(f::SumField, ndt::Int) = steps(f, ndt/min(austrip.(period.((f.a,f.b)))...))

for fun in [:wavelength, :period, :frequency, :wavenumber, :fundamental, :photon_energy]
    @eval begin
        function ($fun)(f::SumField)
            a = ($fun)(f.a)
            b = ($fun)(f.b)
            a != b && error("$(titlecase(string($fun))) differs between SumField composants!")
            a
        end
    end
end

max_frequency(f::SumField) =
    max(max_frequency(f.a), max_frequency(f.b))

continuity(f::SumField) =
    min(continuity(f.a), continuity(f.b))

dimensions(f::SumField) = dimensions(f.a)

phase_shift(f::SumField, ϕ) =
    SumField(phase_shift(f.a, ϕ), phase_shift(f.b, ϕ))

# ** Wrapped fields

"""
    WrappedField

Wrapper around any electric `field`
"""
abstract type WrappedField <: AbstractField end

for fun in [:params, :carrier, :envelope, :polarization,
            :wavelength, :period, :frequency, :max_frequency,
            :wavenumber, :fundamental, :photon_energy,
            :intensity, :amplitude, :duration, :continuity,
            :span, :dimensions]
    @eval $fun(f::WrappedField, args...) = $fun(parent(f), args...)
end
# [:vector_potential, :field_amplitude], should these be explicitly forwarded?

# ** Negated fields

import Base: -

"""
    NegatedField(a)

Represents a field whose [`vector_potential`](@ref) is the negative of
that of `a`.

# Example

```jldoctest
julia> @field(A) do
           I₀ = 1.0
           T = 2.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> B = -A
ElectricFields.NegatedField{ElectricFields.LinearField{ElectricFields.FixedCarrier{Quantity{Float64, 𝐋, Unitful.FreeUnits{(Eₕ⁻¹, ħ, c), 𝐋, nothing}}, Quantity{Float64, 𝐓, Unitful.FreeUnits{(Eₕ⁻¹, ħ), 𝐓, nothing}}, Float64, Int64}, ElectricFields.GaussianEnvelope{Float64}, Float64}}(Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm)

julia> field_amplitude(A, 0.5)
0.008830294277641045

julia> field_amplitude(B, 0.5)
-0.008830294277641045

julia> field_amplitude(A-A, 0.5)
-0.0
```
"""
struct NegatedField{F<:AbstractField} <: WrappedField
    a::F
end

-(a::AbstractField,
  b::AbstractField) = a + NegatedField(b)
-(a::AbstractField) = NegatedField(a)

Base.parent(f::NegatedField) = f.a
vector_potential(f::NegatedField, t) = -vector_potential(parent(f), t)

# ** Delayed fields

"""
    DelayedField(a, t₀)

Represents a delayed copy of `a`, that appears at time `t₀` instead of
at `0`, i.e. `t₀>0` implies the field comes _later_.

# Examples

```jldoctest
julia> @field(A) do
           I₀ = 1.0
           T = 2.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> delay(A, 1u"fs")
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
  – delayed by 41.3414 jiffies = 1.0000 fs

julia> delay(A, 1.0)
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
  – delayed by 1.0000 jiffies = 24.1888 as

julia> delay(A, π*u"rad")
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
  – delayed by 1.0000 jiffies = 24.1888 as
```
"""
struct DelayedField{F<:AbstractField,T} <: WrappedField
    a::F
    t₀::T
end

vector_potential(f::DelayedField, t::Number) = vector_potential(f.a, t-f.t₀)

function show(io::IO, f::DelayedField)
    show(io, f.a)
    printfmt(io, "\n  – delayed by {1:.4f} jiffies = {2:s}",
             f.t₀, au2si_round(f.t₀, u"s"))
end

function span(f::DelayedField)
    s = span(parent(f))
    (s.left+f.t₀)..(s.right+f.t₀)
end

Base.parent(f::DelayedField) = f.a

# *** DONE Delay operators
#     Convention for delayed fields: a field delayed by a /positive/
#     time, comes /later/, i.e. we write \(f(t-\delta t)\).

delay(a::AbstractField, t₀::Union{<:Real,<:Unitful.Time}) = DelayedField(a, austrip(t₀))
delay(a::AbstractField, ϕ::Quantity{<:Real, NoDims}) = delay(a, (ϕ/(2π*u"rad"))*period(a))

delay(a::DelayedField) = a.t₀
delay(a::AbstractField) = 0

export delay

# ** Padded fields

"""
      PaddedField(field, a, b)

Wrapper around any electric `field`, padded with `a` units of time
before, and `b` units of time after the ordinary [`span`](@ref) of the
field.

# Example

```jldoctest
julia> @field(A) do
           I₀ = 1.0
           T = 2.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> B = PaddedField(A, 10.0, 30.0)
Padding before 10.0000 jiffies = 241.8884 as and after 30.0000 jiffies = 725.6653 as of
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> span(A), span(B)
(-6.0..6.0, -16.0..36.0)
```
"""
struct PaddedField{Field<:AbstractField,T} <: WrappedField
    field::Field
    a::T
    b::T

    function PaddedField(field::Field, a::T, b::T) where {Field,T<:Real}
        (a < 0 || b < 0) &&
            throw(ArgumentError("Padding must be non-negative"))
        new{Field,T}(field, a, b)
    end

    PaddedField(field, a::Unitful.Time, b::Unitful.Time) =
        PaddedField(field, austrip(a), austrip(b))
end

Base.parent(f::PaddedField) = f.field

function show(io::IO, f::PaddedField)
    printfmtln(io, "Padding before {1:.4f} jiffies = {2:s} and after {3:.4f} jiffies = {4:s} of",
               f.a, au2si_round(f.a, u"s"),
               f.b, au2si_round(f.b, u"s"))
    show(io, f.field)
end

function vector_potential(f::PaddedField, t::T) where {T<:Number}
    s = span(f.field)
    v = vector_potential(parent(f), clamp(t, endpoints(s)...))
    t ∈ s ? v : zero(v)
end

function span(f::PaddedField)
    a,b = endpoints(span(f.field))
    (a-f.a)..(b+f.b)
end

time_integral(f::PaddedField) =
    time_integral(parent(f))

export PaddedField

# ** Windowed fields

"""
    WindowedField(field, a, b)

Wrapper around any electric `field`, windowed such that it is zero
outside the time interval `a..b`.

# Example

```jldoctest
julia> @field(A) do
           I₀ = 1.0
           T = 2.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> B = WindowedField(A, -3, 5)
Window from -3.0000 jiffies = -72.5665 as to 5.0000 jiffies = 120.9442 as of
Linearly polarized field with
  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
    - A₀ = 0.3183 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm

julia> span(A), span(B)
(-6.0..6.0, -3.0..5.0)

julia> field_amplitude(A, -4)
-0.6395632315635295

julia> field_amplitude(B, -4)
0.0
```
"""
struct WindowedField{Field<:AbstractField,T} <: WrappedField
    field::Field
    a::T
    b::T

    WindowedField(field::Field, a::T, b::T) where {Field,T<:Real} =
        new{Field,T}(field, a, b)

    WindowedField(field, a::Unitful.Time, b::Unitful.Time) =
        WindowedField(field, austrip(a), austrip(b))
end

function show(io::IO, f::WindowedField)
    printfmtln(io, "Window from {1:.4f} jiffies = {2:s} to {3:.4f} jiffies = {4:s} of",
               f.a, au2si_round(f.a, u"s"),
               f.b, au2si_round(f.b, u"s"))
    show(io, f.field)
end

Base.parent(f::WindowedField) = f.field

span(f::WindowedField) = span(parent(f)) ∩ (f.a..f.b)

phase_shift(f::WindowedField, δϕ) =
    WindowedField(phase_shift(parent(f), δϕ), f.a, f.b)

for fun in [:vector_potential, :field_amplitude, :intensity]
    @eval function $fun(f::WindowedField{T}, t) where T
        v = $fun(parent(f), t)
        t < f.a || t > f.b ? zero(v) : v
    end
end

function field_amplitude(f::WindowedField, a, b)
    v = field_amplitude(parent(f), max(a, f.a), min(b, f.b))
    if a ∈ f.a..f.b || b ∈ f.a..f.b
        v
    else
        zero(v)
    end
end

function timeaxis(f::WindowedField, fs)
    t = timeaxis(f.field, fs)
    t[findall(in(f.a..f.b), t)]
end

export WindowedField