function out = get_laserfield(options)
% generate output with pulse parameters, electric field
% Input: tlist [fs], wavelength [nm], intensity [W/cm^2], t0 [fs], phase [pi], width [fs]
% Output: E0 [V/Å], w [rad/fs], T [fs], FWHM [fs], tau [fs], Up [eV], A(t) [V
% fs/Å], E(t) [V/Å], Ephot [eV]
arguments
    options.wavelength (1,1) {mustBeNumeric} = 800
    options.intensity (1,1) {mustBeNumeric} = 1e12
    options.t0 (1,1) {mustBeNumeric} = 0
    options.phase (1,1) {mustBeNumeric} = 0
    options.width (1,1) {mustBeNumeric} = -1
    options.tlist (1,:) double {mustBeNumeric} = 0
end

load_const;

out.E0 = sqrt(2.0*options.intensity*1e4/(const.epsilon_0*const.c))*1e-10;   % electric field [V/Å]
out.w = 2.*pi.*const.c/options.wavelength * 1e-6;                           % angular frequency [1/fs]
out.T = 2*pi/out.w;                                                         % laser period [fs]
if options.width == -1
    out.FWHM = pi/out.w;                                                    % intensity width (FWHM) [fs]
else
    out.FWHM = options.width;
end

if options.tlist == 0
    options.tlist = linspace(options.t0 - 3*out.FWHM, options.t0 + 3*out.FWHM, 1000);
elseif options.tlist(end) < options.t0 + 3*out.FWHM
    warning('time list might be not sufficiently long')
end

out.tau = out.FWHM/(2.*sqrt(log(2)));                                       % width [fs]

out.Up = const.e.^2*(out.E0*1e10).^2./(4*const.m_e*(out.w*1e15).^2*const.e);% ponderomotive potential [eV]
dt = options.tlist(2) - options.tlist(1);
out.t0 = options.t0;
out.phase = options.phase*pi;
out.envelope = exp(-(options.tlist - options.t0).^2/(2.*out.tau.^2));
phi = out.phase + pi/2;
out.A = out.E0/out.w.*cos(out.w.*(options.tlist - options.t0) + phi).*out.envelope;     % vector potential [fs V/Å]
out.E = -gradient(out.A)/dt;

out.Ephot = const.hbar * out.w * 1e15 / const.e;                            % photon energy [eV]

out.tlist = options.tlist;

out.cycles = out.FWHM/pi * out.w;                                           % number of cycles
out.xquiver = const.e*out.E0*1e10/(const.m_e*(out.w*1e15)^2)*1e10;          % quiver aamplitude [Å]

end