RB-SFA.m

(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
(*
This is the RB-SFA package for calculating high-order harmonic generation within the Strong Field Approximation. For the notebook that generated this package file and additional documentaion, see https://github.com/episanty/RB-SFA.
*)


(* ::Input::Initialization:: *)
BeginPackage["RBSFA`"];


(* ::Input::Initialization:: *)
$RBSFAversion::usage="$RBSFAversion prints the current version of the RB-SFA package in use and its timestamp.";
$RBSFAtimestamp::usage="$RBSFAtimestamp prints the timestamp of the current version of the RB-SFA package.";
Begin["`Private`"];
$RBSFAversion:="RB-SFA v2.2.2, "<>$RBSFAtimestamp;
End[];


(* ::Input::Initialization:: *)
RBSFAversion::usage="RBSFAversion[] has been deprecated in favour of $RBSFAversion.";
RBSFAversion::dprc="RBSFAversion[] has been deprecated in favour of $RBSFAversion.";
Begin["`Private`"];
RBSFAversion[]:=(Message[RBSFAversion::dprc];$RBSFAversion);
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];
$RBSFAtimestamp="Tue 6 Oct 2020 18:45:38";
End[];


(* ::Input::Initialization:: *)
$RBSFAdirectory::usage="$RBSFAdirectory is the directory where the current RB-SFA package instance is located.";


(* ::Input::Initialization:: *)
Begin["`Private`"];
With[{softLinkTestString=StringSplit[StringJoin[ReadList["! ls -la "<>StringReplace[$InputFileName,{" "->"\\ "}],String]]," -> "]},
If[Length[softLinkTestString]>1,(*Testing in case $InputFileName is a soft link to the actual directory.*)
$RBSFAdirectory=StringReplace[DirectoryName[softLinkTestString[[2]]],{" "->"\\ "}],
$RBSFAdirectory=StringReplace[DirectoryName[$InputFileName],{" "->"\\ "}];
]];
End[];


(* ::Input::Initialization:: *)
$RBSFAcommit::usage="$RBSFAcommit returns the git commit log at the location of the RB-SFA package if there is one.";
$RBSFAcommit::OS="$RBSFAcommit has only been tested on Linux.";


(* ::Input::Initialization:: *)
Begin["`Private`"];
$RBSFAcommit:=(If[$OperatingSystem!="Unix",Message[$RBSFAcommit::OS]];
StringJoin[Riffle[ReadList["!cd "<>$RBSFAdirectory<>" && git log -1",String],{"\n"}]]);
End[];


(* ::Input::Initialization:: *)
Quiet[Check[
ConstantArray[0,{}];,
Unprotect[ConstantArray];
ConstantArray[Private`x_,{}]:=Private`x;
Protect[ConstantArray];
]];


(* ::Input::Initialization:: *)
Parallelize;
Parallel`Developer`$InitCode=Hold[
Quiet[Check[
ConstantArray[0,{}];,
Unprotect[ConstantArray];
ConstantArray[Private`x_,{}]:=Private`x;
Protect[ConstantArray];
]];
];


(* ::Input::Initialization:: *)
If[
Context[ReIm]=!="System`"&&Attributes[ReIm]=={},
ReIm::usage="\!\(\*RowBox[{\"ReIm\", \"[\", StyleBox[\"z\", \"TI\"], \"]\"}]\) gives the list \!\(\*RowBox[{\"{\", RowBox[{RowBox[{\"Re\", \"[\", StyleBox[\"z\", \"TI\"], \"]\"}], \",\", RowBox[{\"Im\", \"[\", StyleBox[\"z\", \"TI\"], \"]\"}]}], \"}\"}]\) of the number \!\(\*StyleBox[\"z\", \"TI\"]\).";
ReIm[Private`z_]:={Re[Private`z],Im[Private`z]};
SetAttributes[ReIm,Listable];
Protect[ReIm];
]


(* ::Input::Initialization:: *)
AssociationTranspose::usage="AssociationTranspose[association] transposes the given two-level association of associations.";
AssociationTranspose::wrngshp="Input `1` is the wrong shape; it must be an association all of whose Values are valid associations.";
Begin["`Private`"];
AssociationTranspose[association_?(
And@@(AssociationQ/@Join[{#},Values[#]])&
)]:=GroupBy[
Join@@Thread/@Normal//@association,
{First@*Last,First}
][[All,All,1,2,2]];
AssociationTranspose[association__]:="Doesn't display; cf. mm.se/q/29321 for details"/;Message[AssociationTranspose::wrngshp,association]
End[];


(* ::Input::Initialization:: *)
If[
$VersionNumber<10.1,
KeyValueMap::usage="\!\(\*RowBox[{\"KeyValueMap\", \"[\", RowBox[{StyleBox[\"f\",\"TI\"], \",\", RowBox[{\"\[LeftAssociation]\",RowBox[{RowBox[{SubscriptBox[StyleBox[\"key\", \"TI\"], StyleBox[\"1\", \"TR\"]], \"\[Rule]\", SubscriptBox[StyleBox[\"val\", \"TI\"], StyleBox[\"1\", \"TR\"]]}], \",\", RowBox[{SubscriptBox[StyleBox[\"key\", \"TI\"], StyleBox[\"2\", \"TR\"]], \"\[Rule]\", SubscriptBox[StyleBox[\"val\", \"TI\"], StyleBox[\"2\", \"TR\"]]}], \",\", StyleBox[\"\[Ellipsis]\", \"TR\"]}], \"\[RightAssociation]\"}]}], \"]\"}]\) gives the list \!\(\*RowBox[{\"{\", RowBox[{RowBox[{StyleBox[\"f\", \"TI\"], \"[\", RowBox[{SubscriptBox[StyleBox[\"key\", \"TI\"], StyleBox[\"1\", \"TR\"]], \",\", SubscriptBox[StyleBox[\"val\", \"TI\"], StyleBox[\"1\", \"TR\"]]}], \"]\"}], \",\", RowBox[{StyleBox[\"f\", \"TI\"], \"[\", RowBox[{SubscriptBox[StyleBox[\"key\", \"TI\"], StyleBox[\"2\", \"TR\"]], \",\", SubscriptBox[StyleBox[\"val\", \"TI\"], StyleBox[\"2\", \"TR\"]]}], \"]\"}], \",\", StyleBox[\"\[Ellipsis]\", \"TR\"]}], \"}\"}]\). (Note: function backported from v10.1+.)

\!\(\*RowBox[{\"KeyValueMap\", \"[\", StyleBox[\"f\", \"TI\"], \"]\"}]\) represents an operator form of KeyValueMap that can be applied to an expression.";
KeyValueMap::invak="The argument `1` is not a valid association";
]

Begin["`Private`"];
If[
$VersionNumber<10.1,
KeyValueMap[f_,assoc_?AssociationQ]:=Map[Apply[f],Normal[assoc]];
KeyValueMap[f_][assoc_?AssociationQ]:=KeyValueMap[f,assoc];
KeyValueMap[f_,assoc__]:="Doesn't display; cf. mm.se/q/29321 for details"/;Message[KeyValueMap::invak,assoc];
]
End[];


(* ::Input::Initialization:: *)
hydrogenicDTME::usage="hydrogenicDTME[p,\[Kappa]] returns the dipole transition matrix element for a 1s hydrogenic state of ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\).

hydrogenicDTME[p,\[Kappa],{n,l,m}] returns the dipole transition matrix element for an n,l,m hydrogenic state of ground-state ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\).

hydrogenicDTME[p,\[Kappa],n,l,m] returns the dipole transition matrix element for an n,l,m hydrogenic state of ground-state ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\).";
hydrogenicDTMERegularized::usage="hydrogenicDTMERegularized[p,\[Kappa]] returns the dipole transition matrix element for a 1s hydrogenic state of ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\), regularized to remove the denominator of 1/(\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)\!\(\*SuperscriptBox[\()\), \(3\)]\), where the saddle-point solutions are singular.

hydrogenicDTMERegularized[p,\[Kappa],{n,l,m}] returns the dipole transition matrix element for an n,l,m hydrogenic state of ground-state ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\), regularized to remove factors of (\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)) from the denominator.

hydrogenicDTMERegularized[p,\[Kappa],n,l,m] returns the dipole transition matrix element for an n,l,m hydrogenic state of ground-state ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\)=\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\), regularized to remove factors of (\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)) from the denominator.";
Begin["`Private`"];
hydrogenicDTME[p_List,\[Kappa]_]:=(8I)/\[Pi] (Sqrt[2\[Kappa]^5]p)/(Total[p^2]+\[Kappa]^2)^3
hydrogenicDTME[p_?NumberQ,\[Kappa]_]:=(8I)/\[Pi] (Sqrt[2\[Kappa]^5]p)/(p^2+\[Kappa]^2)^3
hydrogenicDTMERegularized[p_List,\[Kappa]_]:=(8I)/\[Pi] (Sqrt[2\[Kappa]^5]p)/1
hydrogenicDTMERegularized[p_?NumberQ,\[Kappa]_]:=(8I)/\[Pi] (Sqrt[2\[Kappa]^5]p)/1
End[];


(* ::Input::Initialization:: *)
gaussianDTME::usage="gaussianDTME[p,\[Kappa]] returns the dipole transition matrix element for a gaussian state of characteristic size 1/\[Kappa].";
Begin["`Private`"];
gaussianDTME[p_List,\[Kappa]_]:=-I (4\[Pi])^(3/4) \[Kappa]^(-7/2) p Exp[-(Total[p^2]/(2\[Kappa]^2))]
gaussianDTME[p_?NumberQ,\[Kappa]_]:=-I (4\[Pi])^(3/4) \[Kappa]^(-7/2) p Exp[-(p^2/(2\[Kappa]^2))]
End[];


(* ::Input::Initialization:: *)
SolidHarmonicS::usage="SolidHarmonicS[l,m,x,y,z] calculates the solid harmonic \!\(\*SubscriptBox[\(S\), \(lm\)]\)(x,y,z)=\!\(\*SuperscriptBox[\(r\), \(l\)]\)\!\(\*SubscriptBox[\(Y\), \(lm\)]\)(x,y,z).

SolidHarmonicS[l,m,{x,y,z}] does the same.";
Begin["`Private`"];
SolidHarmonicS[\[Lambda]_Integer,\[Mu]_Integer,x_,y_,z_]/;\[Lambda]>=Abs[\[Mu]]:=Sqrt[(2 \[Lambda]+1)/(4 \[Pi])] Sqrt[Gamma[\[Lambda]-Abs[\[Mu]]+1]/Gamma[\[Lambda]+Abs[\[Mu]]+1]] 2^-\[Lambda] (-1)^((\[Mu]-Abs[\[Mu]])/2)*
If[Rationalize[\[Mu]]==0,1,(x+Sign[\[Mu]]I y)^Abs[\[Mu]]]*
Sum[
(-1)^(\[Mu]+k) Binomial[\[Lambda],k] Binomial[2 \[Lambda]-2 k,\[Lambda]] Pochhammer[\[Lambda]-Abs[\[Mu]]-2 k+1,Abs[\[Mu]]] *
If[TrueQ[Pochhammer[\[Lambda]-Abs[\[Mu]]-2 k+1,Abs[\[Mu]]]==0],1,
If[Rationalize[k]==0,1,(x^2+y^2+z^2)^k]If[Rationalize[\[Lambda]-Abs[\[Mu]]-2 k]==0,1,z^(\[Lambda]-Abs[\[Mu]]-2 k)]
]
,{k,0,Quotient[\[Lambda],2]}]
SolidHarmonicS[\[Lambda]_Integer,\[Mu]_Integer,{x_,y_,z_}]/;\[Lambda]>=Abs[\[Mu]]:=SolidHarmonicS[\[Lambda],\[Mu],x,y,z]
End[];


(* ::Input::Initialization:: *)
hydrogenic\[CapitalPsi]::usage="hydrogenic\[CapitalPsi][n,l,m,\[Kappa],px,py,pz] calculates the momentum-space wavefunction \[CapitalPsi](p)=\[LeftAngleBracket]p|nlm\[RightAngleBracket] for a hydrogenic atom with ionization potential \!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)/2.

hydrogenic\[CapitalPsi][n,l,m,\[Kappa],{px,py,pz}] calculates the momentum-space wavefunction \[CapitalPsi](p)=\[LeftAngleBracket]p|nlm\[RightAngleBracket] for a hydrogenic atom with ionization potential \!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)/2.";
Begin["`Private`"];
hydrogenic\[CapitalPsi][n_,l_,m_,\[Kappa]\[Kappa]_,ppx_,ppy_,ppz_]:=Block[{\[Kappa],px,py,pz},
hydrogenic\[CapitalPsi][n,l,m,\[Kappa]_,px_,py_,pz_]=Simplify[
-SolidHarmonicS[l,m,px,py,pz] ((-I)^l \[Pi] 2^(2l+4) l!)/(2\[Pi] \[Kappa])^(3/2) Sqrt[(n (n-l-1)!)/(n+l)!] \[Kappa]^(l+4)/(px^2+py^2+pz^2+\[Kappa]^2)^(l+2) GegenbauerC[n-l-1,l+1,(px^2+py^2+pz^2-\[Kappa]^2)/(px^2+py^2+pz^2+\[Kappa]^2)]
];
hydrogenic\[CapitalPsi][n,l,m,\[Kappa]\[Kappa],ppx,ppy,ppz]
];
hydrogenic\[CapitalPsi][n_,l_,m_,\[Kappa]_,{px_,py_,pz_}]:=hydrogenic\[CapitalPsi][n,l,m,\[Kappa],px,py,pz];
End[];


(* ::Input::Initialization:: *)
hydrogenic\[CapitalPsi]Regularized::usage="hydrogenic\[CapitalUpsilon]Regularized[n,l,m,\[Kappa],px,py,pz] calculates the momentum-space wavefunction \[CapitalPsi](p)=\[LeftAngleBracket]p|nlm\[RightAngleBracket] for a hydrogenic atom with ionization potential \!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)/2, multiplied by (\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)\!\(\*SuperscriptBox[\()\), \(n + 1\)]\) to remove any factors of \!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\) in the denominator.

hydrogenic\[CapitalUpsilon]Regularized[n,l,m,\[Kappa],{px,py,pz}] calculates the momentum-space wavefunction \[CapitalPsi](p)=\[LeftAngleBracket]p|nlm\[RightAngleBracket] for a hydrogenic atom with ionization potential \!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)/2, multiplied by (\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)\!\(\*SuperscriptBox[\()\), \(n + 1\)]\) to remove any factors of \!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\) in the denominator.";
Begin["`Private`"];
hydrogenic\[CapitalPsi]Regularized[n_,l_,m_,\[Kappa]\[Kappa]_,ppx_,ppy_,ppz_]:=Block[{\[Kappa],px,py,pz},
hydrogenic\[CapitalPsi]Regularized[n,l,m,\[Kappa]_,px_,py_,pz_]=Simplify[Cancel[
-SolidHarmonicS[l,m,px,py,pz] ((-I)^l \[Pi] 2^(2l+4) l!)/(2\[Pi] \[Kappa])^(3/2) Sqrt[(n (n-l-1)!)/(n+l)!] \[Kappa]^(l+4) (px^2+py^2+pz^2+\[Kappa]^2)^(n-l-1) GegenbauerC[n-l-1,l+1,(px^2+py^2+pz^2-\[Kappa]^2)/(px^2+py^2+pz^2+\[Kappa]^2)]
]];
hydrogenic\[CapitalPsi]Regularized[n,l,m,\[Kappa]\[Kappa],ppx,ppy,ppz]
];
hydrogenic\[CapitalPsi]Regularized[n_,l_,m_,\[Kappa]_,{px_,py_,pz_}]:=hydrogenic\[CapitalPsi]Regularized[n,l,m,\[Kappa],px,py,pz];
End[];


(* ::Input::Initialization:: *)
hydrogenic\[CapitalUpsilon]::usage="hydrogenic\[CapitalUpsilon][n,l,m,\[Kappa],px,py,pz] calculates the Upsilon function \[CapitalUpsilon](p)=(\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SubscriptBox[\(I\), \(p\)]\))\[LeftAngleBracket]p|nlm\[RightAngleBracket] for a hydrogenic atom with ionization potential \!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)/2.

hydrogenic\[CapitalUpsilon][n,l,m,\[Kappa],{px,py,pz}] calculates the Upsilon function \[CapitalUpsilon](p)=(\!\(\*FractionBox[\(1\), \(2\)]\)\!\(\*SuperscriptBox[\(p\), \(2\)]\)+\!\(\*SubscriptBox[\(I\), \(p\)]\))\[LeftAngleBracket]p|nlm\[RightAngleBracket] for a hydrogenic atom with ionization potential \!\(\*SuperscriptBox[\(\[Kappa]\), \(2\)]\)/2.";
Begin["`Private`"];
hydrogenic\[CapitalUpsilon][n_,l_,m_,\[Kappa]_,px_,py_,pz_]:=1/2 (px^2+py^2+pz^2+\[Kappa]^2)hydrogenic\[CapitalPsi][n,l,m,\[Kappa],px,py,pz];
hydrogenic\[CapitalUpsilon][n_,l_,m_,\[Kappa]_,{px_,py_,pz_}]:=hydrogenic\[CapitalUpsilon][n,l,m,\[Kappa],px,py,pz];
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];
hydrogenicDTME[{ppx_,ppy_,ppz_},\[Kappa]\[Kappa]_,n_,l_,m_]:=Block[{\[Kappa],px,py,pz},
hydrogenicDTME[{px_,py_,pz_},\[Kappa]_,n,l,m]=Simplify[Grad[hydrogenic\[CapitalUpsilon][n,l,m,\[Kappa],px,py,pz],{px,py,pz}]];
hydrogenicDTME[{ppx,ppy,ppz},\[Kappa]\[Kappa],n,l,m]
];
hydrogenicDTME[{px_,py_,pz_},\[Kappa]_,{n_,l_,m_}]:=hydrogenicDTME[{px,py,pz},\[Kappa],n,l,m];
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];
hydrogenicDTMERegularized[{px_,py_,pz_},\[Kappa]_,n_,l_,m_]:=(px^2+py^2+pz^2+\[Kappa]^2)^(n+1) hydrogenicDTME[{px,py,pz},\[Kappa],n,l,m];
hydrogenicDTMERegularized[{px_,py_,pz_},\[Kappa]_,{n_,l_,m_}]:=hydrogenicDTMERegularized[{px,py,pz},\[Kappa],n,l,m];
End[];


(* ::Input::Initialization:: *)
flatTopEnvelope::usage="flatTopEnvelope[\[Omega],n,nRamp] returns a Function object representing a flat-top envelope at carrier frequency \[Omega], which lasts for a total of n cycles and is zero afterwards, with a sine-squared ramp on either end lasting for nRamp cycles.";
Begin["`Private`"];
flatTopEnvelope[\[Omega]_,num_,nRamp_]:=Function[t,Piecewise[{{0,t<0},{Sin[(\[Omega] t)/(4nRamp)]^2,0<=t<(2 \[Pi])/\[Omega] nRamp},{1,(2 \[Pi])/\[Omega] nRamp<=t<(2 \[Pi])/\[Omega] (num-nRamp)},{Sin[(\[Omega] ((2 \[Pi])/\[Omega] num-t))/(4nRamp)]^2,(2 \[Pi])/\[Omega] (num-nRamp)<=t<(2 \[Pi])/\[Omega] num},{0,(2 \[Pi])/\[Omega] num<=t}}]]
End[];


(* ::Input::Initialization:: *)
cosPowerFlatTop::usage="cosPowerFlatTop[\[Omega],num,power] returns a Function object representing a smooth flat-top envelope of the form 1-Cos(\[Omega] t/2 num\!\(\*SuperscriptBox[\()\), \(power\)]\)";
Begin["`Private`"];
cosPowerFlatTop[\[Omega]_,num_,power_]:=Function[t,1-Cos[(\[Omega] t)/(2num)]^power]
End[];


(* ::Input::Initialization:: *)
PointsPerCycle::usage="PointsPerCycle is a sampling option which specifies the number of sampling points per cycle to be used in integrations.";
TotalCycles::usage="TotalCycles is a sampling option which specifies the total number of periods to be integrated over.";
CarrierFrequency::usage="CarrierFrequency is a sampling option which specifies the carrier frequency to be used.";
CarrierFrequency::default="Warning: no CarrierFrequency was specified, using \[Omega]=`1` a.u. as the default.";
$DefaultCarrierFrequency::usage="Default CarrierFrequency to use when no explicit option is indicated.";

Protect[PointsPerCycle,TotalCycles,CarrierFrequency];


(* ::Input::Initialization:: *)
standardOptions={PointsPerCycle->90,TotalCycles->1,CarrierFrequency->Automatic,IntegrationPointsPerCycle->Automatic};
$DefaultCarrierFrequency=0.057;


(* ::Input::Initialization:: *)

GetCarrierFrequency::usage="GetCarrierFrequency[OptionValue[CarrierFrequency]] returns OptionValue[CarrierFrequency], unless it's set to Automatic, in which case it returns $DefaultCarrierFrequency and issues a warning.

GetCarrierFrequency[\[Omega]] works for any input.";
Begin["`Private`"];
GetCarrierFrequency[optionvalue_]:=If[
optionvalue===Automatic,
Message[CarrierFrequency::default,$DefaultCarrierFrequency];$DefaultCarrierFrequency,
optionvalue
]
End[];


(* ::Input::Initialization:: *)
harmonicOrderAxis::usage="harmonicOrderAxis[opt\[Rule]value] returns a list of frequencies which can be used as a frequency axis for Fourier transforms, scaled in units of harmonic order, for the provided field duration and sampling options.";
TargetLength::usage="TargetLength is an option for harmonicOrderAxis which specifies the total length required of the resulting list.";
LengthCorrection::usage="LengthCorrection is an option for harmonicOrderAxis which allows for manual correction of the length of the resulting list.";
Protect[LengthCorrection,TargetLength];
Begin["`Private`"];
Options[harmonicOrderAxis]=Join[standardOptions,{TargetLength->Automatic,LengthCorrection->1}];
harmonicOrderAxis::target="Invalid TargetLength option `1`. This must be a positive integer or Automatic.";
harmonicOrderAxis[OptionsPattern[]]:=Module[{num=OptionValue[TotalCycles],npp=OptionValue[PointsPerCycle]},
Piecewise[{
{1/num Range[0.,Round[(npp num+1)/2.]-1+OptionValue[LengthCorrection]],OptionValue[TargetLength]===Automatic},
{Round[(npp num+1)/2.]/num Range[0,OptionValue[TargetLength]-1]/OptionValue[TargetLength],IntegerQ[OptionValue[TargetLength]]&&OptionValue[TargetLength]>=0}
},
Message[harmonicOrderAxis::target,OptionValue["TargetLength"]];Abort[]
]
]
End[];


(* ::Input::Initialization:: *)
frequencyAxis::usage="frequencyAxis[opt\[Rule]value] returns a list of frequencies which can be used as a frequency axis for Fourier transforms, in atomic units of frequency, for the provided field duration and sampling options.";
Begin["`Private`"];
Options[frequencyAxis]=Options[harmonicOrderAxis];
frequencyAxis[options:OptionsPattern[]]:=GetCarrierFrequency[OptionValue[CarrierFrequency]]harmonicOrderAxis[options]
End[];


(* ::Input::Initialization:: *)
timeAxis::usage="timeAxis[opt\[Rule]value] returns a list of times which can be used as a time axis ";
TimeScale::usage="TimeScale is an option for timeAxis which specifies the units the list should use: AtomicUnits by default, or LaserPeriods if required.";
AtomicUnits::usage="AtomicUnits is a value for the option TimeScale of timeAxis which specifies that the times should be in atomic units of time.";
LaserPeriods::usage="LaserPeriods is a value for the option TimeScale of timeAxis which specifies that the times should be in multiples of the carrier laser period.";
Protect[TimeScale,AtomicUnits,LaserPeriods];
Begin["`Private`"];
Options[timeAxis]=standardOptions~Join~{TimeScale->AtomicUnits,PointNumberCorrection->0};
timeAxis::scale="Invalid TimeScale option `1`. Available values are AtomicUnits and LaserPeriods";
timeAxis[OptionsPattern[]]:=Block[{T=2\[Pi]/\[Omega],\[Omega]=GetCarrierFrequency[OptionValue[CarrierFrequency]],num=OptionValue[TotalCycles],npp=OptionValue[PointsPerCycle]},
Piecewise[{
{1,OptionValue[TimeScale]===AtomicUnits},
{1/T,OptionValue[TimeScale]===LaserPeriods}
},
Message[timeAxis::scale,OptionValue[TimeScale]];Abort[]
]*Table[t
,{t,0,num (2\[Pi])/\[Omega],num/(num*npp+OptionValue[PointNumberCorrection]) (2\[Pi])/\[Omega]}
]
]
End[];


(* ::Input::Initialization:: *)
VectorFourier::usage="VectorFourier[array,n] performs a Fourier transform on array only on level n of the array.
VectorFourier[array,{n1,n2,\[Ellipsis],nk}] performs a Fourier transform on array only on the specified levels {n1,n2,\[Ellipsis],nk}.";

Begin["`Private`"];

Options[VectorFourier]=Options[Fourier];
VectorFourier::levelspec="The level specification `1` is not correct.";

VectorFourier[array_,level_?NumericQ,opts:OptionsPattern[]]:=VectorFourier[array,{level},opts]
VectorFourier[array_,levelsToTransform_,opts:OptionsPattern[]]:=Block[{totalLevels,untouchedLevels,permutation,inversePermutation},
totalLevels=Range[Depth[array]-1];
untouchedLevels=Complement[totalLevels,levelsToTransform];
If[Not[SubsetQ[totalLevels,levelsToTransform]],Message[VectorFourier::levelspec,levelsToTransform];Abort[]];
permutation=Join[untouchedLevels,levelsToTransform];
inversePermutation=InversePermutation[permutation];

Transpose[
Map[
Function[Fourier[#,opts]],
Transpose[
array,
inversePermutation
],
{Length[untouchedLevels]}
],
permutation
]
]

End[];


(* ::Input::Initialization:: *)
differentiateDipoleList::usage="differentiateDipoleList[dipoleList,\[Delta]t,DifferentiationOrder\[Rule]n] differentiates dipoleList numerically with step \[Delta]t to order n=0,1,2.";

Begin["`Private`"];

Options[differentiateDipoleList]=Join[{DifferentiationOrder->0},standardOptions];

differentiateDipoleList::diffOrd="Invalid differentiation order `1`.";

differentiateDipoleList[dipoleList_,\[Delta]t_,opts:OptionsPattern[]]:=Block[{},
Piecewise[{
{dipoleList,OptionValue[DifferentiationOrder]==0},
{1/(2\[Delta]t) (Most[Most[dipoleList]]-Rest[Rest[dipoleList]]),OptionValue[DifferentiationOrder]==1},
{1/\[Delta]t^2 (Most[Most[dipoleList]]-2Most[Rest[dipoleList]]+Rest[Rest[dipoleList]]),OptionValue[DifferentiationOrder]==2}},
Message[differentiateDipoleList::diffOrd,OptionValue[DifferentiationOrder]];Abort[]
]
]

End[];


(* ::Input::Initialization:: *)
getSpectralAmplitude::usage="getSpectralAmplitude[DipoleList] returns the spectral amplitude of DipoleList.";

Polarization::usage="Polarization is an option for getSpectrum and getSpectralAmplitude which specifies a polarization vector along which to polarize the dipole list. The default, Polarization\[Rule]False, specifies an unpolarized spectrum.";
ComplexPart::usage="ComplexPart is an option for getSpectrum and getSpectralAmplitude which specifies a function (like Re, Im, or by default #&) which should be applied to the dipole list before the spectrum is taken.";
\[Omega]Power::usage="\[Omega]Power is an option for getSpectrum and getSpectralAmplitude which specifies a power of frequency which should multiply the spectrum.";
DifferentiationOrder::usage="DifferentiationOrder is an option for getSpectrum and getSpectralAmplitude which specifies the order to which the dipole list should be differentiated before the spectrum is taken.";

Protect[Polarization,ComplexPart,\[Omega]Power,DifferentiationOrder];

Begin["`Private`"];

Options[getSpectralAmplitude]=Join[{Polarization->False,ComplexPart->Re,\[Omega]Power->0,DifferentiationOrder->0},standardOptions];

getSpectralAmplitude::\[Omega]Pow="Invalid \[Omega] power `1`.";

getSpectralAmplitude[dipoleList_,OptionsPattern[]]:=Block[
{polarizationVector,preprocessedList,depth,dimensions,
num=OptionValue[TotalCycles],npp=OptionValue[PointsPerCycle],\[Omega],\[Delta]t=(2\[Pi]/\[Omega])/npp
},
polarizationVector=OptionValue[Polarization]/Norm[OptionValue[Polarization]];

preprocessedList=OptionValue[ComplexPart][
differentiateDipoleList[dipoleList,\[Delta]t,DifferentiationOrder->OptionValue[DifferentiationOrder]]
];

If[NumberQ[OptionValue[\[Omega]Power]],Null;,Message[getSpectralAmplitude::\[Omega]Pow,OptionValue[\[Omega]Power]];Abort[]  ];
If[OptionValue[\[Omega]Power]!=0,\[Omega]=GetCarrierFrequency[OptionValue[CarrierFrequency]],\[Omega]=1];
(*If \[Omega]Power\[Equal]0 the value of \[Omega] doesn't matter and there's no sense in risking printing error messages that would result from a missing CarrierFrequency option.*)

Times[
Sqrt[num] Table[(\[Omega]/num k)^OptionValue[\[Omega]Power],{k,1,Round[Length[preprocessedList]/2]}],
Map[
If[
OptionValue[Polarization]===False,
#&,(*unpolarized spectrum*)
Dot[#,polarizationVector]&
],
VectorFourier[preprocessedList,{1}]
][[1;;Round[Length[preprocessedList]/2]]]
]

]

End[];


(* ::Input::Initialization:: *)
getSpectrum::usage="getSpectrum[DipoleList] returns the power spectrum of DipoleList.";


Begin["`Private`"];
Options[getSpectrum]=Options[getSpectralAmplitude];

getSpectrum[dipoleList_,opts:OptionsPattern[]]:=Map[
Norm[#]^2&,
getSpectralAmplitude[dipoleList,opts]
]

End[];


(* ::Input::Initialization:: *)
spectrumPlotter::usage="spectrumPlotter[spectrum] plots the given spectrum with an appropriate axis in a \!\(\*SubscriptBox[\(log\), \(10\)]\) scale.";
FrequencyAxis::usage="FrequencyAxis is an option for spectrumPlotter which specifies the axis to use.";
Protect[FrequencyAxis];
Begin["`Private`"];
Options[spectrumPlotter]=Join[{FrequencyAxis->"HarmonicOrder"},Options[harmonicOrderAxis],Options[ListLinePlot]];
spectrumPlotter[spectrum_,options:OptionsPattern[]]:=ListPlot[
{Which[
OptionValue[FrequencyAxis]==="HarmonicOrder",
harmonicOrderAxis["TargetLength"->Length[spectrum],Sequence@@FilterRules[{options}~Join~Options[spectrumPlotter],Options[harmonicOrderAxis]]],
OptionValue[FrequencyAxis]==="Frequency",
frequencyAxis["TargetLength"->Length[spectrum],Sequence@@FilterRules[{options}~Join~Options[spectrumPlotter],Options[harmonicOrderAxis]]],
True,Range[Length[spectrum]]
],
Log[10,spectrum]
}\[Transpose]
,Sequence@@FilterRules[{options},Options[ListLinePlot]]
,Joined->True
,PlotRange->Full
,PlotStyle->Thick
,Frame->True
,Axes->False
,ImageSize->800
]
End[];


(* ::Input::Initialization:: *)
biColorSpectrum::usage="biColorSpectrum[DipoleList] produces a two-colour spectrum of DipoleList, separating the two circular polarizations.";
Begin["`Private`"];
Options[biColorSpectrum]=Join[{PlotRange->All},Options[Show],Options[spectrumPlotter],DeleteCases[Options[getSpectrum],Polarization->False]];
biColorSpectrum[dipoleList_,options:OptionsPattern[]]:=Show[{
spectrumPlotter[
getSpectrum[dipoleList,Polarization->{1,+I},Sequence@@FilterRules[{options},Options[getSpectrum]]],
PlotStyle->Red,Sequence@@FilterRules[{options},Options[spectrumPlotter]]],
spectrumPlotter[
getSpectrum[dipoleList,Polarization->{1,-I},Sequence@@FilterRules[{options},Options[getSpectrum]]],
PlotStyle->Blue,Sequence@@FilterRules[{options},Options[spectrumPlotter]]]
}
,PlotRange->OptionValue[PlotRange]
,Sequence@@FilterRules[{options},Options[Show]]
]
End[];


(* ::Input::Initialization:: *)
SineSquaredGate::usage="SineSquaredGate[nGateRamp] specifies an integration gate with a sine-squared ramp, such that SineSquaredGate[nGateRamp][\[Omega]t,nGate] has nGate flat periods and nGateRamp ramp periods.";
LinearRampGate::usage="LinearRampGate[nGateRamp] specifies an integration gate with a linear ramp, such that SineSquaredGate[nGateRamp][\[Omega]t,nGate] has nGate flat periods and nGateRamp ramp periods.";
Begin["`Private`"];
SineSquaredGate[nGateRamp_][\[Omega]\[Tau]_,nGate_]:=Piecewise[{{1,\[Omega]\[Tau]<=2\[Pi] (nGate-nGateRamp)},{Sin[(2\[Pi] nGate-\[Omega]\[Tau])/(4nGateRamp)]^2,2\[Pi] (nGate-nGateRamp)<\[Omega]\[Tau]<=2\[Pi] nGate},{0,nGate<\[Omega]\[Tau]}}]
LinearRampGate[nGateRamp_][\[Omega]\[Tau]_,nGate_]:=Piecewise[{{1,\[Omega]\[Tau]<=2\[Pi] (nGate-nGateRamp)},{-((\[Omega]\[Tau]-2\[Pi] (nGate+nGateRamp))/(2\[Pi] nGateRamp)),2\[Pi] (nGate-nGateRamp)<\[Omega]\[Tau]<=2\[Pi] nGate},{0,nGate<\[Omega]\[Tau]}}]
End[];


(* ::Input::Initialization:: *)
getIonizationPotential::usage="getIonizationPotential[Target] returns the ionization potential of an atomic target, e.g. \"Hydrogen\", in atomic units.\[IndentingNewLine]
getIonizationPotential[Target,q] returns the ionization potential of the q-th ion of the specified Target, in atomic units.\[IndentingNewLine]
getIonizationPotential[{Target,q}] returns the ionization potential of the q-th ion of the specified Target, in atomic units.";
Begin["`Private`"];
getIonizationPotential[Target_,Charge_:0]:=UnitConvert[ElementData[Target,"IonizationEnergies"][[Charge+1]]/(Quantity[1,"AvogadroConstant"]Quantity[1,"Hartrees"])]
getIonizationPotential[{Target_,Charge_:0}]:=getIonizationPotential[Target,Charge]
End[];


(* ::Input::Initialization:: *)
makeDipoleList::usage="makeDipoleList[VectorPotential\[Rule]A] calculates the dipole response to the vector potential A.";

VectorPotential::usage="VectorPotential is an option for makeDipole list which specifies the field's vector potential. Usage should be VectorPotential\[Rule]A, where A[t]//.pars must yield a list of numbers for numeric t and parameters indicated by FieldParameters\[Rule]pars.";
VectorPotentialGradient::usage="VectorPotentialGradient is an option for makeDipole list which specifies the gradient of the field's vector potential. Usage should be VectorPotentialGradient\[Rule]GA, where GA[t]//.pars must yield a square matrix of the same dimension as the vector potential for numeric t and parameters indicated by FieldParameters\[Rule]pars. The indices must be such that GA[t]\[LeftDoubleBracket]i,j\[RightDoubleBracket] returns \!\(\*SubscriptBox[\(\[PartialD]\), \(i\)]\)\!\(\*SubscriptBox[\(A\), \(j\)]\)[t].";
ElectricField::usage="ElectricField is an option for makeDipole list which specifies an electric field to use in the ionization matrix element, in case the time derivative of the vector potential is not desired. Usage should be ElectricField\[Rule]F, where F[t]//.pars must yield a list of numbers for numeric t and parameters indicated by FieldParameters\[Rule]pars.";
FieldParameters::usage="FieldParameters is an option for makeDipole list which provides the parameters, as a list of replacement rules, to use in evaluating the vector potential. If the vector potential is provided as VectorPotential\[Rule]A, and the parameters as FieldParameters\[Rule]pars, then A[t]//.pars must yield a list of numbers if given a numeric argument t.";
Preintegrals::usage="Preintegrals is an option for makeDipole list which specifies whether the preintegrals of the vector potential should be \"Analytic\" or \"Numeric\".";
ReportingFunction::usage="ReportingFunction is an option for makeDipole list which specifies a function used to report the results, either internally (by the default, Identity) or to an external file.";
Gate::usage="Gate is an option for makeDipole list which specifies the integration gate to use. Usage as Gate\[Rule]g, nGate\[Rule]n will gate the integral at time \[Omega]t/\[Omega] by g[\[Omega]t,n]. The default is Gate\[Rule]SineSquaredGate[1/2].";
nGate::usage="nGate is an option for makeDipole list which specifies the total number of cycles in the integration gate.";
IonizationPotential::usage="IonizationPotential is an option for makeDipoleList which specifies the ionization potential \!\(\*SubscriptBox[\(I\), \(p\)]\) of the target.";
Target::usage="Target is an option for makeDipoleList which specifies chemical species producing the HHG emission, pulling the ionization potential from the Wolfram ElementData curated data set.";
DipoleTransitionMatrixElement::usage="DipoleTransitionMatrixElement is an option for makeDipoleList which secifies a function f to use as the dipole transition matrix element, or a pair of functions {\!\(\*SubscriptBox[\(f\), \(ion\)]\),\!\(\*SubscriptBox[\(f\), \(rec\)]\)} to be used separately for the ionization and recombination dipoels, to be used in the form f[p,\[Kappa]]=f[p,\!\(\*SqrtBox[\(2 \*SubscriptBox[\(I\), \(p\)]\)]\)].";
\[Epsilon]Correction::usage="\[Epsilon]Correction is an option for makeDipoleList which specifies the regularization correction \[Epsilon], i.e. as used in the factor \!\(\*FractionBox[\(1\), SuperscriptBox[\((t - tt + \[ImaginaryI]\[Epsilon])\), \(3/2\)]]\).";
PointNumberCorrection::usage="PointNumberCorrection is an option for makeDipoleList and timeAxis which specifies an extra number of points to be integrated over, which is useful to prevent Indeterminate errors when a Piecewise envelope is being differentiated at the boundaries.";
IntegrationPointsPerCycle::usage="IntegrationPointsPerCycle is an option for makeDipoleList which controls the number of points per cycle to use for the integration. Set to Automatic, to follow PointsPerCycle, or to an integer.";
RunInParallel::usage="RunInParallel is an option for makeDipoleList which controls whether each RB-SFA instance is parallelized. It accepts False as the (Automatic) option, True, to parallelize each instance, or a pair of functions {TableCommand, SumCommand} to use for the iteration and summing, which could be e.g. {Inactive[ParallelTable], Inactive[Sum]}.";
Simplifier::usage="Simplifier is an option for makeDipoleList which specifies a function to use to simplify the intermediate and final analytical results.";
CheckNumericFields::usage="CheckNumericFields is an option for makeDipoleList which specifies whether to check for numeric values of A[t] and GA[t] for numeric t.";
QuadraticActionTerms::usage="QuadraticActionTerms is an option for makeDipoleList which specifies whether to use quadratic terms in \[Del]\!\(\*SuperscriptBox[\(A\), \(2\)]\) in the action.";

Protect[VectorPotential,VectorPotentialGradient,ElectricField,FieldParameters,Preintegrals,ReportingFunction,Gate,nGate,IonizationPotential,Target,\[Epsilon]Correction,PointNumberCorrection,DipoleTransitionMatrixElement,IntegrationPointsPerCycle,RunInParallel,Simplifier,CheckNumericFields,QuadraticActionTerms];



Begin["`Private`"];
Options[makeDipoleList]=standardOptions~Join~{
VectorPotential->Automatic,FieldParameters->{},VectorPotentialGradient->None,ElectricField->Automatic,
Preintegrals->"Analytic",ReportingFunction->Identity,
Gate->SineSquaredGate[1/2],nGate->3/2,\[Epsilon]Correction->0.1,
IonizationPotential->0.5,Target->Automatic,DipoleTransitionMatrixElement->hydrogenicDTME,
PointNumberCorrection->0,Verbose->0,CheckNumericFields->True,
RunInParallel->Automatic,
Simplifier->Identity,QuadraticActionTerms->True
};
makeDipoleList::gate="The integration gate g provided as Gate\[Rule]`1` is incorrect. Its usage as g[`2`,`3`] returns `4` and should return a number.";
makeDipoleList::pot="The vector potential A provided as VectorPotential\[Rule]`1` is incorrect or is missing FieldParameters. Its usage as A[`2`] returns `3` and should return a list of numbers.";
makeDipoleList::efield="The electric field f provided as ElectricField\[Rule]`1` is incorrect or is missing FieldParameters. Its usage as F[`2`] returns `3` and should return a list of numbers. Alternatively, use ElectricField\[Rule]Automatic.";
makeDipoleList::gradpot="The vector potential GA provided as VectorPotentialGradient\[Rule]`1` is incorrect or is missing FieldParameters. Its usage as GA[`2`] returns `3` and should return a square matrix of numbers. Alternatively, use VectorPotentialGradient\[Rule]None.";
makeDipoleList::preint="Wrong Preintegrals option `1`. Valid options are \"Analytic\" and \"Numeric\".";
makeDipoleList::runpar="Wrong RunInParallel option `1`.";
makeDipoleList::carrfreq="Non-numeric option CarrierFrequency `1`.";




makeDipoleList[OptionsPattern[]]:=Block[
{
num=OptionValue[TotalCycles],npp=OptionValue[PointsPerCycle],\[Omega],
dipoleRec,dipoleIon,\[Kappa],
A,F,GA,pi,ps,S,
gate,tGate,setPreintegral,
tInit,tFinal,\[Delta]t,\[Delta]tint,\[Epsilon]=OptionValue[\[Epsilon]Correction],
AInt,A2Int,GAInt,GAdotAInt,AdotGAInt,GAIntInt,
PScorrectionInt,constCorrectionInt,GAIntdotGAIntInt,QuadMatrix,q,
simplifier,prefactor,integrand,dipoleList,
TableCommand,SumCommand
},

A[t_]=OptionValue[VectorPotential][t]//.OptionValue[FieldParameters];
If[
OptionValue[ElectricField]===Automatic,F[t_]=-D[A[t],t];,
F[t_]=OptionValue[ElectricField][t]//.OptionValue[FieldParameters];
];
GA[t_]=If[
TrueQ[OptionValue[VectorPotentialGradient]==None],        Table[0,{Length[A[tInit]]},{Length[A[tInit]]}],
OptionValue[VectorPotentialGradient][t]//.OptionValue[FieldParameters]
];

\[Omega]=GetCarrierFrequency[OptionValue[CarrierFrequency]];
If[!NumberQ[\[Omega]]&&TrueQ[OptionValue[CheckNumericFields]],Message[makeDipoleList::carrfreq,\[Omega]];Abort[]];
tInit=0;
tFinal=(2\[Pi])/\[Omega] num;
(*looping timestep*)
\[Delta]t=(tFinal-tInit)/(num*npp+OptionValue[PointNumberCorrection]);
(*integration timestep*)
\[Delta]tint=If[OptionValue[IntegrationPointsPerCycle]===Automatic,\[Delta]t,(tFinal-tInit)/(num*OptionValue[IntegrationPointsPerCycle]+OptionValue[PointNumberCorrection])];
tGate=OptionValue[nGate] (2\[Pi])/\[Omega];
(*Check potential and potential gradient for correctness.*)
(*To do: change logic conditions to constructions on VectorQ[#,NumberQ]& and MatrixQ.*)
If[TrueQ[OptionValue[CheckNumericFields]],
With[{\[Omega]tRandom=RandomReal[{\[Omega] tInit,\[Omega] tFinal}]},
If[!And@@(NumberQ/@A[\[Omega]tRandom/\[Omega]]),Message[makeDipoleList::pot,OptionValue[VectorPotential],\[Omega]tRandom,A[\[Omega]tRandom]];Abort[]];
If[!And@@(NumberQ/@Flatten[GA[\[Omega]tRandom/\[Omega]]]),Message[makeDipoleList::gradpot,OptionValue[VectorPotentialGradient],\[Omega]tRandom,GA[\[Omega]tRandom]];Abort[]];
If[!And@@(NumberQ/@F[\[Omega]tRandom/\[Omega]]),Message[makeDipoleList::efield,OptionValue[ElectricField],\[Omega]tRandom,F[\[Omega]tRandom]];Abort[]];
]];

gate[\[Omega]\[Tau]_]:=OptionValue[Gate][\[Omega]\[Tau],OptionValue[nGate]];
With[{\[Omega]tRandom=RandomReal[{\[Omega] tInit,\[Omega] tFinal}]},
If[!TrueQ[NumberQ[gate[\[Omega]tRandom]]],
Message[makeDipoleList::gate,OptionValue[Gate],\[Omega]tRandom,OptionValue[nGate],gate[\[Omega]tRandom]];Abort[]]
];

(*Target setup*)
Which[
OptionValue[Target]===Automatic,\[Kappa]=Sqrt[2OptionValue[IonizationPotential]],
True,\[Kappa]=Sqrt[2getIonizationPotential[OptionValue[Target]]]
];
With[{dim=Length[A[RandomReal[{\[Omega] tInit,\[Omega] tFinal}]]]},
(*Explicit conjugation of the recombination matrix element to keep the integrand analytic.*)
Which[
Head[OptionValue[DipoleTransitionMatrixElement]]===List,
dipoleIon[{p1_,p2_,p3_}[[1;;dim]],\[Kappa]\[Kappa]_]=First[OptionValue[DipoleTransitionMatrixElement]][{p1,p2,p3}[[1;;dim]],\[Kappa]\[Kappa]];
dipoleRec[{p1_,p2_,p3_}[[1;;dim]],\[Kappa]\[Kappa]_]=Assuming[{{p1,p2,p3,\[Kappa]\[Kappa]}\[Element]Reals},Simplify[
Conjugate[Last[OptionValue[DipoleTransitionMatrixElement]][{p1,p2,p3}[[1;;dim]],\[Kappa]\[Kappa]]]
]];
,True,
dipoleIon[{p1_,p2_,p3_}[[1;;dim]],\[Kappa]\[Kappa]_]=OptionValue[DipoleTransitionMatrixElement][{p1,p2,p3}[[1;;dim]],\[Kappa]\[Kappa]];
dipoleRec[{p1_,p2_,p3_}[[1;;dim]],\[Kappa]\[Kappa]_]=Assuming[{{p1,p2,p3,\[Kappa]\[Kappa]}\[Element]Reals},Simplify[
Conjugate[OptionValue[DipoleTransitionMatrixElement][{p1,p2,p3}[[1;;dim]],\[Kappa]\[Kappa]]]
]];
];
];


simplifier=OptionValue[Simplifier];
q=Boole[TrueQ[OptionValue[QuadraticActionTerms]]];

setPreintegral[integralVariable_,preintegrand_,dimensions_,integrateWithoutGradient_,parametric_]:=Which[
OptionValue[VectorPotentialGradient]=!=None||TrueQ[integrateWithoutGradient],(*Vector potential gradient specified, or integral variable does not depend on it, so integrate*)
Which[
OptionValue[Preintegrals]=="Analytic",
integralVariable[t_,tt_]=simplifier[((#/.{\[Tau]->t})-(#/.{\[Tau]->tt}))&[Integrate[preintegrand[\[Tau],tt],\[Tau]]]];

,OptionValue[Preintegrals]=="Numeric",
Which[
TrueQ[Not[parametric]],
Block[{innerVariable},
integralVariable[t_,tt_]=(innerVariable[t]-innerVariable[tt]/.First[
NDSolve[{innerVariable'[\[Tau]]==preintegrand[\[Tau]],innerVariable[tInit]==ConstantArray[0,dimensions]},innerVariable,{\[Tau],tInit,tFinal},MaxStepSize->0.25/\[Omega]]
])
];
,True,
Block[{matrixpreintegrand,innerVariable,\[Tau]pre},
matrixpreintegrand[indices_,t_?NumericQ,tt_?NumericQ]:=preintegrand[t,tt][[##&@@indices]];
integralVariable[t_,tt_]=Array[(
innerVariable[##][t-tt,tt]/.First@NDSolve[{
D[innerVariable[##][\[Tau]pre,tt],\[Tau]pre]==Piecewise[{{matrixpreintegrand[{##},tt+\[Tau]pre,tt],tt+\[Tau]pre<=tFinal}},0],
innerVariable[##][0,tt]==0
},innerVariable[##]
,{\[Tau]pre,0,tFinal-tInit},{tt,tInit,tFinal}
,MaxStepSize->0.25/\[Omega]
]
)&,dimensions];
]
];
];
,OptionValue[VectorPotentialGradient]===None,(*Vector potential gradient has not been specified, and integral variable depends on it, so return appropriate zero matrix*)
integralVariable[t_]=ConstantArray[0,dimensions];
integralVariable[t_,tt_]=ConstantArray[0,dimensions];
];

Apply[setPreintegral,({
 {AInt, A[#1]&, {Length[A[tInit]]}, True, False},
 {A2Int, A[#1].A[#1]&, {}, True, False},
 {GAInt, GA[#1]&, {Length[A[tInit]],Length[A[tInit]]}, False, False},
 {GAdotAInt, GA[#1].A[#1]&, {Length[A[tInit]]}, False, False},
 {AdotGAInt, A[#1].GA[#1]&, {Length[A[tInit]]}, False, False},
 {GAIntInt, GAInt[#1,#2]&, {Length[A[tInit]],Length[A[tInit]]}, False, True},
 {PScorrectionInt, GAdotAInt[#1,#2]+A[#1].GAInt[#1,#2]-q GAInt[#1,#2]\[Transpose].GAdotAInt[#1,#2]&, {Length[A[tInit]]}, False, True},
 {GAIntdotGAIntInt, q GAInt[#1,#2]\[Transpose].GAInt[#1,#2]&, {Length[A[tInit]],Length[A[tInit]]}, False, True},
 {constCorrectionInt, (A[#1]-q/2 GAdotAInt[#1,#2]).GAdotAInt[#1,#2]&, {}, False, True}
}),{1}];
(*{\!\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]\(A\((\[Tau])\)\[DifferentialD]\[Tau]\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]\(A
\*SuperscriptBox[\((\[Tau])\), \(2\)]\[DifferentialD]\[Tau]\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]\(\[Del]A\((\[Tau])\)\[DifferentialD]\[Tau]\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]\(\[Del]A\((\[Tau])\)\[CenterDot]A\((\[Tau])\)\[DifferentialD]\[Tau]\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]\(A\((\[Tau])\)\[CenterDot]\[Del]A\((\[Tau])\)\[DifferentialD]\[Tau]\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(t\)]\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]\[Del]A\((\[Tau]')\)\[DifferentialD]\[Tau]'\[DifferentialD]\[Tau]\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(t\)]\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]
\*SubscriptBox[\(\[PartialD]\), \(j\)]
\*SubscriptBox[\(A\), \(k\)]\((\[Tau]')\)
\*SubscriptBox[\(A\), \(k\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\)\)+Subscript[A, k](\[Tau])\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]\(
\*SubscriptBox[\(\[PartialD]\), \(k\)]
\*SubscriptBox[\(A\), \(j\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\)\)-\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]\(
\*SubscriptBox[\(\[PartialD]\), \(i\)]
\*SubscriptBox[\(A\), \(j\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]
\*SubscriptBox[\(\[PartialD]\), \(i\)]
\*SubscriptBox[\(A\), \(k\)]\((\[Tau]')\)
\*SubscriptBox[\(A\), \(k\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\[DifferentialD]\[Tau]\)\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(t\)]\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]
\*SubscriptBox[\(\[PartialD]\), \(i\)]
\*SubscriptBox[\(A\), \(j\)]\((\[Tau]')\)
\*SubscriptBox[\(A\), \(j\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\(
\*SubsuperscriptBox[\(\[Integral]\), 
SubscriptBox[\(t\), \(0\)], \(t\)]
\*SubscriptBox[\(\[PartialD]\), \(i\)]
\*SubscriptBox[\(A\), \(k\)]\((\[Tau]')\)
\*SubscriptBox[\(A\), \(k\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\[DifferentialD]\[Tau]\)\)\),\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(t\)]\(\((
\*SubscriptBox[\(A\), \(k\)]\((\[Tau])\) - 
\*FractionBox[\(1\), \(2\)]\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]
\*SubscriptBox[\(\[PartialD]\), \(k\)]
\*SubscriptBox[\(A\), \(i\)]\((\[Tau]')\)
\*SubscriptBox[\(A\), \(i\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\))\)\[CenterDot]\(
\*SubsuperscriptBox[\(\[Integral]\), \(t'\), \(\[Tau]\)]
\*SubscriptBox[\(\[PartialD]\), \(k\)]
\*SubscriptBox[\(A\), \(j\)]\((\[Tau]')\)
\*SubscriptBox[\(A\), \(j\)]\((\[Tau]')\)\[DifferentialD]\[Tau]'\[DifferentialD]\[Tau]\)\)\)};*)


(*Displaced momentum*)
pi[p_,t_,tt_]:=p+A[t]-GAInt[t,tt].p-GAdotAInt[t,tt];

(*Quadratic coefficient in nondipole action*)
QuadMatrix[t_,tt_]:=(GAIntInt[t,tt]+GAIntInt[t,tt]\[Transpose])/2-1/2 GAIntdotGAIntInt[t,tt];

(*Stationary momentum and action*)
ps[t_,tt_]:=ps[t,tt]=-(1/(t-tt-I \[Epsilon]))Inverse[IdentityMatrix[Length[A[tInit]]]-1/(t-tt-I \[Epsilon]) 2QuadMatrix[t,tt]].(AInt[t,tt]-PScorrectionInt[t,tt]);
S[t_,tt_]:=simplifier[
1/2 (Total[ps[t,tt]^2]+\[Kappa]^2)(t-tt)+ps[t,tt].AInt[t,tt]+1/2 A2Int[t,tt]-(
ps[t,tt].QuadMatrix[t,tt].ps[t,tt]+ps[t,tt].PScorrectionInt[t,tt]+constCorrectionInt[t,tt]
)
];

prefactor[t_,\[Tau]_]:=I ((2\[Pi])/(\[Epsilon]+I \[Tau]))^(3/2) dipoleRec[pi[ps[t,t-\[Tau]],t,t-\[Tau]],\[Kappa]]*dipoleIon[pi[ps[t,t-\[Tau]],t-\[Tau],t-\[Tau]],\[Kappa]].F[t-\[Tau]];
integrand[t_,\[Tau]_]:=prefactor[t,\[Tau]]Exp[-I S[t,t-\[Tau]]]gate[\[Omega] \[Tau]];


(*Debugging constructs. Verbose\[Rule]1 prints information about the internal functions. Verbose\[Rule]2 returns all the relevant internal functions and stops. Verbose\[Rule]3 for quantum-orbit constructs.*)
Which[
OptionValue[Verbose]==1,Information/@{A,GA,ps,pi,S,AInt,A2Int,GAInt,GAdotAInt,AdotGAInt,GAIntInt,PScorrectionInt,constCorrectionInt,GAIntdotGAIntInt},
OptionValue[Verbose]==2,Return[With[{t=Symbol["t"],tt=Symbol["tt"],\[Tau]=Symbol["\[Tau]"],p={Symbol["p1"],Symbol["p2"],Symbol["p3"]}[[1;;Length[A[\[Omega] tInit]]]]},
{A[t],GA[t],ps[t,tt],pi[p,t,tt],S[t,tt],AInt[t,tt],A2Int[t,tt],GAInt[t,tt],GAdotAInt[t,tt],AdotGAInt[t,tt],GAIntInt[t,tt],PScorrectionInt[t,tt],constCorrectionInt[t,tt],GAIntdotGAIntInt[t,tt],QuadMatrix[t,tt],integrand[t,\[Tau]]}]],
OptionValue[Verbose]==3,
Return[{
Function[Evaluate[prefactor[#1,#1-#2]]],Function[Evaluate[S[#1,#2]]]
}]
];

(*Single-run parallelization*)
Which[
OptionValue[RunInParallel]===Automatic||OptionValue[RunInParallel]===False, TableCommand=Table;SumCommand=Sum;,
OptionValue[RunInParallel]===True,TableCommand=ParallelTable;SumCommand=Sum;,
True,TableCommand=OptionValue[RunInParallel][[1]];SumCommand=OptionValue[RunInParallel][[2]];
];

(*Numerical integration loop*)
dipoleList=Table[
OptionValue[ReportingFunction][
\[Delta]tint Sum[(
integrand[t,\[Tau]]
),{\[Tau],0,If[OptionValue[Preintegrals]=="Analytic",tGate,Min[t-tInit,tGate]],\[Delta]tint}]
]
,{t,tInit,tFinal,\[Delta]t}
];
dipoleList

]
End[];


(* ::Input::Initialization:: *)
FindComplexRoots::usage="FindComplexRoots[e1==e2, {z, zmin, zmax}] attempts to find complex roots of the equation e1==e2 in the complex rectangle with corners zmin and zmax.

FindComplexRoots[{e1==e2, e3==e4, \[Ellipsis]}, {z1, z1min, z1max}, {z2, z2min, z2max}, \[Ellipsis]] attempts to find complex roots of the given system of equations in the multidimensional complex rectangle with corners z1min, z1max, z2min, z2max, \[Ellipsis].";
Seeds::usage="Seeds is an option for FindComplexRoots which determines how many initial seeds are used to attempt to find roots of the given equation.";
SeedGenerator::usage="SeedGenerator is an option for FindComplexRoots which determines the function   used to generate the seeds for the internal FindRoot call. Its value can be RandomComplex, RandomNiederreiterComplexes, RandomSobolComplexes, DeterministicComplexGrid, or any function f such that f[{zmin, zmax}, n] returns n complex numbers in the rectancle with corners zmin and zmax.";

Options[FindComplexRoots] = Join[Options[FindRoot], {Seeds -> 50, SeedGenerator -> RandomComplex, Tolerance -> Automatic, Verbose -> False}];
SyntaxInformation[FindComplexRoots] = {"ArgumentsPattern" -> {_, {_, _, _}, OptionsPattern[]},    "LocalVariables" -> {"Table", {2, \[Infinity]}}};
FindComplexRoots::seeds = "Value of option Seeds -> `1` is not a positive integer.";
FindComplexRoots::tol = "Value of option Tolerance -> `1` is not Automatic or a number in [0,\[Infinity]).";
$MessageGroups=Join[$MessageGroups,{"FindComplexRoots":>{FindRoot::lstol}}];

Protect[Seeds];
Protect[SeedGenerator];


(* ::Input::Initialization:: *)
SetTolerances::usage="SetTolerances[tolerance,length] produces a list of the given length with the specified tolerance, which may be a number or a list of numbers.\n
SetTolerances[tolerance,length,workingPrecision] allows a fallback to a specified workingPrecision in case the given tolerance fails to be numeric.";
Begin["`Private`"];
SetTolerances[tolerance_,length_,workingPrecision_:$MachinePrecision]:=Which[
ListQ[tolerance],tolerance,
True,ConstantArray[
Which[
NumberQ[tolerance],tolerance,
True,10^If[NumberQ[workingPrecision], 2-workingPrecision,2-$MachinePrecision]
]
,length]
]
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];
FindComplexRoots[equations_List,domainSpecifiers__, ops : OptionsPattern[]] := Block[{seeds,tolerances},
If[! IntegerQ[Rationalize[OptionValue[Seeds]]] || OptionValue[Seeds]<=0,Message[FindComplexRoots::seeds, OptionValue[Seeds]]];If[! (OptionValue[Tolerance] === Automatic || OptionValue[Tolerance]>=0),Message[FindComplexRoots::tol, OptionValue[Seeds]]];

seeds=OptionValue[SeedGenerator][{domainSpecifiers}[[All,{2,3}]],OptionValue[Seeds]];
tolerances=SetTolerances[OptionValue[Tolerance],Length[{domainSpecifiers}],OptionValue[WorkingPrecision]];

If[OptionValue[Verbose],Hold[], Hold[FindRoot::lstol]] /. {
Hold[messageSequence___] :> Quiet[
DeleteDuplicates[
Select[
Check[
FindRoot[
equations
,Evaluate[Sequence@@Table[{{domainSpecifiers}[[j,1]],#[[j]]},{j,Length[{domainSpecifiers}]}]]
,Evaluate[Sequence @@ FilterRules[{ops}, Options[FindRoot]]]
],
##&[]
]&/@seeds,
Function[
repList,
ReplaceAll[
Evaluate[And@@Table[
And[
Re[{domainSpecifiers}[[j,2]]]<=Re[{domainSpecifiers}[[j,1]]]<=Re[{domainSpecifiers}[[j,3]]],
Im[{domainSpecifiers}[[j,2]]]<=Im[{domainSpecifiers}[[j,1]]]<=Im[{domainSpecifiers}[[j,3]]]
]
,{j,Length[{domainSpecifiers}]}]]
,repList]
]
],
Function[{repList1,repList2},
And@@Table[
Abs[({domainSpecifiers}[[j,1]]/.repList1)-({domainSpecifiers}[[j,1]]/.repList2)]<tolerances[[j]]
,{j,Length[{domainSpecifiers}]}]
]
]
, {messageSequence}]}
]
FindComplexRoots[e1_==e2_,{z_,zmin_,zmax_},ops:OptionsPattern[]]:=FindComplexRoots[{e1==e2},{z,zmin,zmax},ops]
End[];


(* ::Input::Initialization:: *)
RandomSobolComplexes::usage="RandomSobolComplexes[{zmin, zmax}, n] generates a low-discrepancy Sobol sequence of n quasirandom complex numbers in the rectangle with corners zmin and zmax.

RandomSobolComplexes[{{z1min,z1max},{z2min,z2max},\[Ellipsis]},n] generates a low-discrepancy Sobol sequence of n quasirandom complex numbers in the multi-dimensional rectangle with corners {z1min,z1max},{z2min,z2max},\[Ellipsis].";


(* ::Input::Initialization:: *)
Begin["`Private`"];
RandomSobolComplexes[pairsList__, number_] :=Map[
Function[randomsList,
pairsList[[All,1]]+Complex@@@Times[
ReIm[pairsList[[All,2]]-pairsList[[All,1]]],
randomsList
]
],
BlockRandom[
SeedRandom[Method->{"MKL",Method->{"Sobol", "Dimension" -> 2Length[pairsList]}}];
SeedRandom[];
RandomReal[{0, 1}, {number,Length[pairsList],2}]
]
]
RandomSobolComplexes[{zmin_?NumericQ,zmax_?NumericQ},number_]:=RandomSobolComplexes[{{zmin,zmax}},number][[All,1]]
End[];


(* ::Input::Initialization:: *)
RandomNiederreiterComplexes::usage="RandomNiederreiterComplexes[{zmin, zmax}, n] generates a low-discrepancy Niederreiter sequence of n quasirandom complex numbers in the rectangle with corners zmin and zmax.

RandomNiederreiterComplexes[{{z1min,z1max},{z2min,z2max},\[Ellipsis]},n] generates a low-discrepancy Niederreiter sequence of n quasirandom complex numbers in the multi-dimensional rectangle with corners {z1min,z1max},{z2min,z2max},\[Ellipsis].";


(* ::Input::Initialization:: *)
Begin["`Private`"];
RandomNiederreiterComplexes[pairsList__, number_] :=Map[
Function[randomsList,
pairsList[[All,1]]+Complex@@@Times[
ReIm[pairsList[[All,2]]-pairsList[[All,1]]],
randomsList
]
],
BlockRandom[
SeedRandom[Method->{"MKL",Method->{"Niederreiter", "Dimension" -> 2Length[pairsList]}}];
SeedRandom[];
RandomReal[{0, 1}, {number,Length[pairsList],2}]
]
]
RandomNiederreiterComplexes[{zmin_?NumericQ,zmax_?NumericQ},number_]:=RandomNiederreiterComplexes[{{zmin,zmax}},number][[All,1]]
End[];


(* ::Input::Initialization:: *)
DeterministicComplexGrid::usage="DeterministicComplexGrid[{zmin, zmax}, n] generates a grid of about n equally spaced complex numbers in the rectangle with corners zmin and zmax.

DeterministicComplexGrid[{{z1min,z1max},{z2min,z2max},\[Ellipsis]},n] generates a regular grid of about n equally spaced complex numbers in the multi-dimensional rectangle with corners {z1min,z1max},{z2min,z2max},\[Ellipsis].";


(* ::Input::Initialization:: *)
Begin["`Private`"];
DeterministicComplexGrid[pairsList_,number_]:=Block[{sep,separationsList,gridPointBasis,k},
sep=NestWhile[0.99#&,Min[Flatten[ReIm[pairsList[[All,2]]-pairsList[[All,1]]]]],Times@@(Floor[Flatten[ReIm[pairsList[[All,2]]-pairsList[[All,1]]]],0.99#]/(0.99#))<=number&];
separationsList=Round[Floor[Flatten[ReIm[pairsList[[All,2]]-pairsList[[All,1]]]],sep]/sep];
gridPointBasis=MapThread[
Function[{l,n},Range[l[[1]],l[[2]],(l[[2]]-l[[1]])/(n+1)][[2;;-2]]],
{Flatten[Transpose[ReIm[pairsList],{1,3,2}],1],separationsList}
];
Flatten[Table[
Table[k[2j-1]+I k[2j],{j,1,Length[pairsList]}],
Evaluate[Sequence@@Table[{k[j],gridPointBasis[[j]]},{j,1,2Length[pairsList]}]]
],Evaluate[Range[1,2Length[pairsList]]]]
]
DeterministicComplexGrid[{zmin_?NumericQ,zmax_?NumericQ},number_]:=DeterministicComplexGrid[{{zmin,zmax}},number][[All,1]]
End[];


(* ::Input::Initialization:: *)
Begin["`Private`"];
Unprotect[RandomComplex];
RandomComplex[{range1_List,moreRanges___},number_]:=Transpose[RandomComplex[#,number]&/@{range1,moreRanges}]
Protect[RandomComplex];
End[];


(* ::Input::Initialization:: *)
Parallelize;
If[Head[Parallel`Developer`$InitCode]=!=Hold,
Parallel`Developer`$InitCode=Hold[]
];
Parallel`Developer`$InitCode=Join[
Parallel`Developer`$InitCode,
Hold[
Unprotect[RandomComplex];
RandomComplex[{Private`range1_List,Private`moreRanges___},Private`number_]:=Transpose[RandomComplex[#,Private`number]&/@{Private`range1,Private`moreRanges}];
Protect[RandomComplex];
]
];


(* ::Input::Initialization:: *)
ConstrainedDerivative::usage="ConstrainedDerivative[n][f][t,tt] calculates the nth derivative of f[t,tt] with respect to t under the constraint that Derivative[0,1][f][t,tt]\[Congruent]0.";

Begin["`Private`"];
ConstrainedDerivative[n_][F_][te_,tte_]:=Block[{f,tts,t,tt},
ConstrainedDerivative[n][f_][t_,tt_]=Nest[
Function[
Simplify[
D[#/.{tt->tts[t]},t]/.{Derivative[0,1][f][t,tts[t]]->0,tts'[t]->-(Derivative[1,1][f][t,tts[t]]/Derivative[0,2][f][t,tts[t]])}
]/.{tts[t]->tt}
]
,f[t,tt],n];
ConstrainedDerivative[n][F][te,tte]
]
End[];


(* ::Input::Initialization:: *)
GetSaddlePoints::usage="GetSaddlePoints[\[CapitalOmega],S,{tmin,tmax},{\[Tau]min,\[Tau]max}] finds a list of solutions {t,\[Tau]} of the HHG temporal saddle-point equations at harmonic energy \[CapitalOmega] for action S, in the range {tmin, tmax} of recombination time and {\[Tau]min, \[Tau]max} of excursion time, where both ranges should be the lower-left and upper-right corners of rectangles in the complex plane.

GetSaddlePoints[\[CapitalOmega]Range,S,{tmin,tmax},{\[Tau]min,\[Tau]max}] finds solutions of the HHG temporal saddle-point equations for a range of harmonic energies \[CapitalOmega]Range, and returns an Association with each harmonic energy \[CapitalOmega] indexing a list of saddle-point solution pairs {t,\[Tau]}.

GetSaddlePoints[\[CapitalOmega]spec,S,{{{\!\(\*SubscriptBox[\(tmin\), \(1\)]\),\!\(\*SubscriptBox[\(tmax\), \(1\)]\)},{\!\(\*SubscriptBox[\(\[Tau]min\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]max\), \(1\)]\)}},{{\!\(\*SubscriptBox[\(tmin\), \(2\)]\),\!\(\*SubscriptBox[\(tmax\), \(2\)]\)},{\!\(\*SubscriptBox[\(\[Tau]min\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]max\), \(2\)]\)}},\[Ellipsis]}] uses multiple time domains and combines the solutions.

GetSaddlePoints[\[CapitalOmega]spec,S,{{urange,vrange},\[Ellipsis]},IndependentVariables\[Rule]{u,v}] uses the explicit independent variables u and v to solve the equations and over the given ranges, where u and v can be any of \"RecombinationTime\", \"IonizationTime\" and \"ExcursionTime\", or their shorthands \"t\", \"tt\" and \"\[Tau]\" resp.";
SortingFunction::usage="SortingFunction is an option of GetSaddlePoints which sets a function f, to be used as f[t,\[Tau],S,\[CapitalOmega]], to be used to sort the solutions, or a list of such functions.";
SelectionFunction::usage="SelectionFunction is an option of GetSaddlePoints that sets a function f, to be used as f[t,\[Tau],S,\[CapitalOmega]], such that roots are only kept if f returns True.";
IndependentVariables::usage="IndependentVariables is an option for GetSaddlePoints that specifies the two independent variables, out of \"RecombinationTime\", \"IonizationTime\" and \"ExcursionTime\" (or their shorthands \"t\", \"tt\" and \"\[Tau]\", respectively), to be used in solving the saddle-point equations, and which range over the given regions.";

FiniteDifference::usage="FiniteDifference is a value for the option Jacobian of FindRoot, FindComplexRoots, GetSaddlePoints, and related functions, which specifies that the Jacobian at each step should be evaluated using numerical finite difference procedures.";

GetSaddlePoints::error="Errors encountered for harmonic energy \[CapitalOmega]=`1`.";

Begin["`Private`"];
Options[GetSaddlePoints]=Join[{SortingFunction->(#2&),SelectionFunction->(True&),IndependentVariables->{"RecombinationTime","ExcursionTime"}},Options[FindComplexRoots]];
Protect[SortingFunction,SelectionFunction,IndependentVariables,FiniteDifference];

GetSaddlePoints[\[CapitalOmega]spec_,S_,{tmin_,tmax_},{\[Tau]min_,\[Tau]max_},options:OptionsPattern[]]:=GetSaddlePoints[\[CapitalOmega]spec,S,{{{tmin,tmax},{\[Tau]min,\[Tau]max}}},options]

GetSaddlePoints[\[CapitalOmega]_,S_,timeRanges_,options:OptionsPattern[]]:=Block[{equations,roots,t=Symbol["t"],tt=Symbol["tt"],\[Tau]=Symbol["\[Tau]"],indVars,depVar,depVarRule,tolerances},
tolerances=SetTolerances[OptionValue[Tolerance],2,OptionValue[WorkingPrecision]];
indVars=OptionValue[IndependentVariables]/.{"RecombinationTime"->"t","ExcursionTime"->"\[Tau]","IonizationTime"->"tt"};
depVar=First[DeleteCases[{"t","\[Tau]","tt"},Alternatives@@indVars]];
depVarRule=depVar/.{"tt"->{tt->t-\[Tau]},"t"->{t->tt+\[Tau]},"\[Tau]"->{\[Tau]->t-tt}};
equations={D[S[t,tt],t]==\[CapitalOmega],D[S[t,tt],tt]==0}/.depVarRule;

SortBy[
DeleteDuplicates[
Flatten[Table[
Select[
Check[
roots=({t,\[Tau]}/.depVarRule)/.(FindComplexRoots[
equations
,Evaluate[Sequence[{Symbol[indVars[[1]]],range[[1,1]],range[[1,2]]},{Symbol[indVars[[2]]],range[[2,1]],range[[2,2]]}]]
,Evaluate[Sequence@@FilterRules[{options},Options[FindComplexRoots]]]
,SeedGenerator->RandomSobolComplexes
,Seeds->50
]/.{{}->(({t,\[Tau]}/.depVarRule)->{})})(*to deal with empty results*)
,Message[GetSaddlePoints::error,\[CapitalOmega]];roots
]
,Function[timesPair,OptionValue[SelectionFunction][timesPair[[1]],timesPair[[2]],S,\[CapitalOmega]]]
]
,{range,timeRanges}],1]
,Function[{timesPair1,timesPair2},    And@@Thread[Abs[timesPair1-timesPair2]<tolerances]     ]
]
,If[
ListQ[OptionValue[SortingFunction]],
Table[Function[timesPair,f[timesPair[[1]],timesPair[[2]],S,\[CapitalOmega]]],{f,OptionValue[SortingFunction]}],
Function[timesPair,OptionValue[SortingFunction][timesPair[[1]],timesPair[[2]],S,\[CapitalOmega]]]
]
]
]
GetSaddlePoints[\[CapitalOmega]Range_List,S_,timeRanges_,options:OptionsPattern[]]:=Association[ParallelTable[
\[CapitalOmega]->GetSaddlePoints[\[CapitalOmega],S,timeRanges,options]
,{\[CapitalOmega],Sort[\[CapitalOmega]Range]}]]

End[];


(* ::Input::Initialization:: *)
GetSaddlesFromSeeds::usage="GetSaddlesFromSeeds[{{\!\(\*SubscriptBox[\(t\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)},{\!\(\*SubscriptBox[\(t\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(2\)]\)},\[Ellipsis]},\[CapitalOmega],S] finds a list of solutions {t,\[Tau]} of the HHG temporal saddle-point equations at harmonic energy \[CapitalOmega] for action S, using the given {\!\(\*SubscriptBox[\(t\), \(i\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(i\)]\)} as seeds for the process.

GetSaddlesFromSeeds[\[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]{{\!\(\*SubscriptBox[\(t\), \(11\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(11\)]\)},{\!\(\*SubscriptBox[\(t\), \(12\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(12\)]\)},\[Ellipsis]},\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\)\[RightArrow]{{\!\(\*SubscriptBox[\(t\), \(21\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(21\)]\)},{\!\(\*SubscriptBox[\(t\), \(22\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(22\)]\)},\[Ellipsis]},\[Ellipsis]\[RightAssociation],\[CapitalOmega],S] finds solutions of the HHG temporal saddle-point equations, using the seeds list from the \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(i\)]\) that's closest to \[CapitalOmega], or as specified by the value of KeyChooserFunction.

GetSaddlesFromSeeds[seeds,{\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\),\[Ellipsis]},S] iterates over the given set of harmonic energies.";

SeedsChooserFunction::usage="SeedsChooserFunction is an option for GetSaddlesFromSeeds that specifies a function f (set by default to Nearest) that, when used as f[{\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\),\[Ellipsis]},\[CapitalOmega]], should return the indices {\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(i\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(j\)]\),\[Ellipsis]} corresponding to the seed sets {{{\!\(\*SubscriptBox[\(t\), \(i1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(i1\)]\)},\[Ellipsis]},{{\!\(\*SubscriptBox[\(t\), \(j1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(j1\)]\)},\[Ellipsis]}} to be used to solve the HHG saddle-point equations.";
RecalculateRoots::usage="RecalculateRoots is an option for GetSaddlesFromSeeds that specifies whether to re-solve the saddle-point equations if the given harmonic energy \[CapitalOmega] is among the set of keys of the given seeds association. The default is False, which is appropriate for S being the same action used to find the seeds, in which case setting RecalculateRoots\[RightArrow]True will produce multiple FindRoot errors. If using a different action than used to find the seeds, set to True.";

GetSaddlesFromSeeds::error="Errors encountered for harmonic energy \[CapitalOmega]=`1`.";
GetSaddlesFromSeeds::norecalc="Skipping re-calculation of roots at harmonic energy `1` since it is already in the key set of the given seeds association. To run the calculation for this case set RecalculateRoots to True.";

Begin["`Private`"];
Options[GetSaddlesFromSeeds]=Join[{RecalculateRoots->False,SeedsChooserFunction->Nearest},Options[GetSaddlePoints]];
Protect[SeedsChooserFunction,RecalculateRoots];


GetSaddlesFromSeeds[seedsSpec_,\[CapitalOmega]Range_List,S_,options:OptionsPattern[]]:=Association[ParallelTable[
\[CapitalOmega]->GetSaddlesFromSeeds[seedsSpec,\[CapitalOmega],S,options]
,{\[CapitalOmega],Sort[\[CapitalOmega]Range]}]]

GetSaddlesFromSeeds[seedsAssociation_Association,\[CapitalOmega]_,S_,options:OptionsPattern[]]:=With[{keys=OptionValue[SeedsChooserFunction][Keys[seedsAssociation],\[CapitalOmega]]},
If[MemberQ[keys,\[CapitalOmega]]&&TrueQ[!OptionValue[RecalculateRoots]],Message[GetSaddlesFromSeeds::norecalc,\[CapitalOmega]];Return[seedsAssociation[\[CapitalOmega]]]];
GetSaddlesFromSeeds[Flatten[Values[seedsAssociation[[Key/@keys]]],1],\[CapitalOmega],S,options]
]


GetSaddlesFromSeeds[seedsList_List,\[CapitalOmega]_?NumberQ,S_,options:OptionsPattern[]]:=Block[
{equations,roots,t=Symbol["t"],tt=Symbol["tt"],\[Tau]=Symbol["\[Tau]"],indVars,depVar,depVarRule,fullSeedVars,tolerances},
tolerances=SetTolerances[OptionValue[Tolerance],2,OptionValue[WorkingPrecision]];
indVars=OptionValue[IndependentVariables]/.{"RecombinationTime"->"t","ExcursionTime"->"\[Tau]","IonizationTime"->"tt"};
depVar=First[DeleteCases[{"t","\[Tau]","tt"},Alternatives@@indVars]];
depVarRule=depVar/.{"tt"->{tt->t-\[Tau]},"t"->{t->tt+\[Tau]},"\[Tau]"->{\[Tau]->t-tt}};
fullSeedVars[seed_]:=<|"t"->seed[[1]],"\[Tau]"->seed[[2]],"tt"->seed[[1]]-seed[[2]]|>;
equations={D[S[t,tt],t]==\[CapitalOmega],D[S[t,tt],tt]==0}/.depVarRule;

SortBy[
DeleteDuplicates[
Select[
Table[
Check[
roots=({t,\[Tau]}/.depVarRule)/.(
FindRoot[
equations
,{Symbol[#],fullSeedVars[seed][[#]]}&/@indVars
,Evaluate[Sequence@@FilterRules[{options},Options[FindRoot]]]
]
/.{{}->(({t,\[Tau]}/.depVarRule)->{})})
,Message[GetSaddlesFromSeeds::error,\[CapitalOmega]];roots
]
,{seed,seedsList}]
,Function[timesPair,OptionValue[SelectionFunction][timesPair[[1]],timesPair[[2]],S,\[CapitalOmega]]]
]
,Function[{timesPair1,timesPair2},    And@@Thread[Abs[timesPair1-timesPair2]<tolerances]     ]
]
,If[
ListQ[OptionValue[SortingFunction]],
Table[Function[timesPair,f[timesPair[[1]],timesPair[[2]],S,\[CapitalOmega]]],{f,OptionValue[SortingFunction]}],
Function[timesPair,OptionValue[SortingFunction][timesPair[[1]],timesPair[[2]],S,\[CapitalOmega]]]
]
]
]
End[];


(* ::Input::Initialization:: *)
GetDoubleSaddlePoints::usage="GetDoubleSaddlePoints[S,{tmin,tmax},{\[Tau]min,\[Tau]max},{\[CapitalOmega]min,\[CapitalOmega]max}] finds a list of double solutions {t,\[Tau],\[CapitalOmega]} of the HHG temporal saddle-point equations, using a complex-valued \[CapitalOmega] in the range {\[CapitalOmega]min,\[CapitalOmega]max}, for action S, in the range {tmin, tmax} of recombination time and {\[Tau]min, \[Tau]max} of excursion time, where the ranges should indicate the lower-left and upper-right corners of rectangles in the complex plane.

GetDoubleSaddlePoints[S,{{{\!\(\*SubscriptBox[\(tmin\), \(1\)]\),\!\(\*SubscriptBox[\(tmax\), \(1\)]\)},{\!\(\*SubscriptBox[\(\[Tau]min\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]max\), \(1\)]\)},{\!\(\*SubscriptBox[\(\[CapitalOmega]min\), \(1\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]max\), \(1\)]\)}},{{\!\(\*SubscriptBox[\(tmin\), \(2\)]\),\!\(\*SubscriptBox[\(tmax\), \(2\)]\)},{\!\(\*SubscriptBox[\(\[Tau]min\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]max\), \(2\)]\)},{\!\(\*SubscriptBox[\(\[CapitalOmega]min\), \(2\)]\),\!\(\*SubscriptBox[\(\[CapitalOmega]max\), \(2\)]\)}},\[Ellipsis]}] uses multiple variable ranges and combines the solutions.

GetDoubleSaddlePoints[S,{{urange,vrange,\[CapitalOmega]range},\[Ellipsis]},IndependentVariables\[Rule]{u,v}] uses the explicit independent temporal variables u and v to solve the equations and over the given ranges, where u and v can be any of \"RecombinationTime\", \"IonizationTime\" and \"ExcursionTime\", or their shorthands \"t\", \"tt\" and \"\[Tau]\" resp.";

GetDoubleSaddlePoints::error="Errors encountered at tag `1`.";

ErrorReportingTag::usage="";
Protect[ErrorReportingTag];


Begin["`Private`"];
Options[GetDoubleSaddlePoints]=Join[{SortingFunction->({#4&,#2&,#1&}),SelectionFunction->(True&),ErrorReportingTag->None,IndependentVariables->{"RecombinationTime","ExcursionTime"}},Options[FindComplexRoots]];

GetDoubleSaddlePoints[S_,{tmin_,tmax_},{\[Tau]min_,\[Tau]max_},{\[CapitalOmega]min_,\[CapitalOmega]max_},options:OptionsPattern[]]:=GetDoubleSaddlePoints[S,{{{tmin,tmax},{\[Tau]min,\[Tau]max},{\[CapitalOmega]min,\[CapitalOmega]max}}},options]

GetDoubleSaddlePoints[S_,variableRanges_,options:OptionsPattern[]]:=Block[{equations,roots,t=Symbol["t"],tt=Symbol["tt"],\[Tau]=Symbol["\[Tau]"],\[CapitalOmega],indVars,depVar,depVarRule,tolerances},
tolerances=SetTolerances[OptionValue[Tolerance],3,OptionValue[WorkingPrecision]];
indVars=OptionValue[IndependentVariables]/.{"RecombinationTime"->"t","ExcursionTime"->"\[Tau]","IonizationTime"->"tt"};
depVar=First[DeleteCases[{"t","\[Tau]","tt"},Alternatives@@indVars]];
depVarRule=depVar/.{"tt"->{tt->t-\[Tau]},"t"->{t->tt+\[Tau]},"\[Tau]"->{\[Tau]->t-tt}};
equations={
D[S[t,tt],t]==\[CapitalOmega],
D[S[t,tt],tt]==0,
D[S[t,tt],{t,2}] D[S[t,tt],{tt,2}]-D[S[t,tt],t,tt]^2==0
}/.depVarRule;

SortBy[
DeleteDuplicates[
Flatten[Table[
Select[
Check[
roots=({t,\[Tau],\[CapitalOmega]}/.depVarRule)/.(FindComplexRoots[
equations
,Evaluate[Sequence[
{Symbol[indVars[[1]]],range[[1,1]],range[[1,2]]},
{Symbol[indVars[[2]]],range[[2,1]],range[[2,2]]},
{\[CapitalOmega],range[[3,1]],range[[3,2]]}
]]
,Evaluate[Sequence@@FilterRules[{options},Options[FindComplexRoots]]]
,SeedGenerator->RandomSobolComplexes
,Seeds->50
]/.{{}->(({t,\[Tau]}/.depVarRule)->{})})(*to deal with empty results*)
,
If[OptionValue[ErrorReportingTag]!=None,Message[GetDoubleSaddlePoints::error,OptionValue[ErrorReportingTag]]];
roots
]
,Function[variableSet,OptionValue[SelectionFunction][variableSet[[1]],variableSet[[2]],S]]
]
,{range,variableRanges}],1]
,Function[{variableSet1,variableSet2},    And@@Thread[Abs[variableSet1-variableSet2]<tolerances]     ]
]
,If[
ListQ[OptionValue[SortingFunction]],
Table[Function[variableSet,f[variableSet[[1]],variableSet[[2]],S,variableSet[[3]]]],{f,OptionValue[SortingFunction]}],
Function[variableSet,OptionValue[SortingFunction][variableSet[[1]],variableSet[[2]],S,variableSet[[3]]]]
]
]
]
End[];


(* ::Input::Initialization:: *)
GetCutoffSaddlePoints::usage="GetCutoffSaddlePoints[S,{tmin,tmax},{\[Tau]min,\[Tau]max}] finds a list of solutions {t,\[Tau],\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)S,\!\(\*SubsuperscriptBox[\(d\), \(t\), \(3\)]\)S} of the HHG cutoff saddle-point equations {\!\(\*SubscriptBox[\(\[PartialD]\), \(t'\)]\)S=0,\!\(\*SubsuperscriptBox[\(d\), \(t\), \(2\)]\)S=0} for action S, in the range {tmin, tmax} of recombination time and {\[Tau]min, \[Tau]max} of excursion time, where both ranges should be the lower-left and upper-right corners of rectangles in the complex plane.

GetCutoffSaddlePoints[S,{{{\!\(\*SubscriptBox[\(tmin\), \(1\)]\),\!\(\*SubscriptBox[\(tmax\), \(1\)]\)},{\!\(\*SubscriptBox[\(\[Tau]min\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]max\), \(1\)]\)}},{{\!\(\*SubscriptBox[\(tmin\), \(2\)]\),\!\(\*SubscriptBox[\(tmax\), \(2\)]\)},{\!\(\*SubscriptBox[\(\[Tau]min\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]max\), \(2\)]\)}},\[Ellipsis]}] uses multiple time domains and combines the solutions.

GetCutoffSaddlePoints[S,{{urange,vrange},\[Ellipsis]},IndependentVariables\[Rule]{u,v}] uses the explicit independent variables u and v to solve the equations and over the given ranges, where u and v can be any of \"RecombinationTime\", \"IonizationTime\" and \"ExcursionTime\", or their shorthands \"t\", \"tt\" and \"\[Tau]\" resp.";

GetCutoffSaddlePoints::error="Errors encountered at tag `1`.";

Begin["`Private`"];
Options[GetCutoffSaddlePoints]=Join[{SortingFunction->(#2&),SelectionFunction->(True&),ErrorReportingTag->None,IndependentVariables->{"RecombinationTime","ExcursionTime"}},Options[FindComplexRoots]];

GetCutoffSaddlePoints[S_,{tmin_,tmax_},{\[Tau]min_,\[Tau]max_},options:OptionsPattern[]]:=GetCutoffSaddlePoints[S,{{{tmin,tmax},{\[Tau]min,\[Tau]max}}},options]

GetCutoffSaddlePoints[S_,timeRanges_,options:OptionsPattern[]]:=Block[{equations,roots,t=Symbol["t"],tt=Symbol["tt"],\[Tau]=Symbol["\[Tau]"],indVars,depVar,depVarRule,tolerances,d1S,d3S,d2Sdtt2},
tolerances=SetTolerances[OptionValue[Tolerance],2,OptionValue[WorkingPrecision]];
indVars=OptionValue[IndependentVariables]/.{"RecombinationTime"->"t","ExcursionTime"->"\[Tau]","IonizationTime"->"tt"};
depVar=First[DeleteCases[{"t","\[Tau]","tt"},Alternatives@@indVars]];
depVarRule=depVar/.{"tt"->{tt->t-\[Tau]},"t"->{t->tt+\[Tau]},"\[Tau]"->{\[Tau]->t-tt}};
equations={
D[S[t,tt],tt]==0,
D[S[t,tt],{t,2}]D[S[t,tt],{tt,2}]-D[S[t,tt],t,tt]^2==0
}/.depVarRule;
d1S[t_,tt_]=ConstrainedDerivative[1][S][t,tt];
d3S[t_,tt_]=ConstrainedDerivative[3][S][t,tt];
d2Sdtt2[t_,tt_]=Derivative[0,2][S][t,tt];

Map[
Association[Thread[{"t","\[Tau]","S","\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)S","\!\(\*SubsuperscriptBox[\(d\), \(t\), \(3\)]\)S","\!\(\*SubsuperscriptBox[\(\[PartialD]\), \(tt\), \(2\)]\)S"}->#]]&,
SortBy[
DeleteDuplicates[
Flatten[Table[
Select[
Check[
roots=({t,\[Tau],S[t,t-\[Tau]],d1S[t,t-\[Tau]],d3S[t,t-\[Tau]],d2Sdtt2[t,t-\[Tau]]}/.depVarRule)/.(FindComplexRoots[
equations
,Evaluate[Sequence[{Symbol[indVars[[1]]],range[[1,1]],range[[1,2]]},{Symbol[indVars[[2]]],range[[2,1]],range[[2,2]]}]]
,Evaluate[Sequence@@FilterRules[{options},Options[FindComplexRoots]]]
,SeedGenerator->RandomSobolComplexes
,Seeds->50
]/.{{}->(({t,\[Tau]}/.depVarRule)->{})})(*to deal with empty results*)
,
If[OptionValue[ErrorReportingTag]!=None,Message[GetCutoffSaddlePoints::error,OptionValue[ErrorReportingTag]]];
roots
]
,Function[timesPair,OptionValue[SelectionFunction][timesPair[[1]],timesPair[[2]],S]]
]
,{range,timeRanges}],1]
,Function[{timesPair1,timesPair2},    And@@Thread[Abs[timesPair1-timesPair2][[{1,2}]]<tolerances]     ]
]
,If[
ListQ[OptionValue[SortingFunction]],
Table[Function[timesPair,f[timesPair[[1]],timesPair[[2]],S]],{f,OptionValue[SortingFunction]}],
Function[timesPair,OptionValue[SortingFunction][timesPair[[1]],timesPair[[2]],S]]
]
]
]
]
End[];


(* ::Input::Initialization:: *)
ClassifyQuantumOrbits::usage="ClassifyQuantumOrbits[saddlePoints,f] sorts an indexed set of saddle points of the form \[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]{{\!\(\*SubscriptBox[\(t\), \(11\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(11\)]\)},{\!\(\*SubscriptBox[\(t\), \(12\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(12\)]\)},\[Ellipsis]}\[Ellipsis]\[RightAssociation] using a function f, which should turn f[t,\[Tau],\[CapitalOmega]] into an appropriate label, and returns an association of the form \[LeftAssociation]\!\(\*SubscriptBox[\(label\), \(1\)]\)\[RightArrow]\[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]\[LeftAssociation]1\[RightArrow]{t,\[Tau]},2\[RightArrow]{t,\[Tau]},\[Ellipsis]\[RightAssociation],\[Ellipsis]\[RightAssociation],\[Ellipsis]\[RightAssociation].

ClassifyQuantumOrbits[saddlePoints,f,sortFunction] uses the function sortFunction to sort the sets of saddle points {{\!\(\*SubscriptBox[\(t\), \(11\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(11\)]\)},{\!\(\*SubscriptBox[\(t\), \(12\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(12\)]\)},\[Ellipsis]} for each label and harmonic energy.

ClassifyQuantumOrbits[saddlePoints,f,sortFunction,DiscardedLabels\[RightArrow]{\!\(\*SubscriptBox[\(label\), \(1\)]\),\!\(\*SubscriptBox[\(label\), \(2\)]\),\[Ellipsis]}] specifies a list of labels to discard from the final output.";
DiscardedLabels::usage="DiscardedLabels is an option for ClassifyQuantumOrbits which specifies a list of labels to discard from the final output.";

Begin["`Private`"];

Options[ClassifyQuantumOrbits]={DiscardedLabels->{}};
Protect[DiscardedLabels];

ClassifyQuantumOrbits[saddlePointList_,classifierFunction_,sortingFunction_:Sort,OptionsPattern[]]:=Map[
Composition[
Association,
MapIndexed[#2[[1]]->#1&],
sortingFunction
],
Delete[
AssociationTranspose[
MapIndexed[
GroupBy[classifierFunction@@#&][Flatten/@Transpose[{#1,ConstantArray[#2[[{1},1]],Length[#1]]}]]&
,saddlePointList][[All,All,All,{1,2}]]
]
,List/@OptionValue[DiscardedLabels]]
,{2}]

End[];


(* ::Input::Initialization:: *)
ClearAll[ReperiodSaddles]
ReperiodSaddles::usage="ReperiodSaddles[{{\!\(\*SubscriptBox[\(t\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)},{\!\(\*SubscriptBox[\(t\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(2\)]\)},\[Ellipsis]},f] readjusts the assigned cycle of the saddle points {\!\(\*SubscriptBox[\(t\), \(i\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(i\)]\)}, returning the list {{\!\(\*SubscriptBox[\(t\), \(1\)]\)+f[\!\(\*SubscriptBox[\(t\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)],\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)},\[Ellipsis]}.

ReperiodSaddles[\[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]{{\!\(\*SubscriptBox[\(t\), \(11\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(11\)]\)},\[Ellipsis]},\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\)\[Rule]\[Ellipsis]\[RightAssociation],f] reperiods saddle-point pairs in a harmonic-energy-indexed association.

ReperiodSaddles[\[LeftAssociation]\!\(\*SubscriptBox[\(label\), \(l\)]\)\[Rule]\[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]{{\!\(\*SubscriptBox[\(t\), \(11\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(11\)]\)},\[Ellipsis]},\[Ellipsis]\[RightAssociation],\[Ellipsis]\[RightAssociation],f] reperiods saddle-point pairs of a classified set of saddle points.";

Begin["`Private`"];

ReperiodSaddles[pair_/;Depth[pair]==2,f_]:={pair[[1]]+f[pair[[1]],pair[[2]]],pair[[2]]}
ReperiodSaddles[association_,f_]:=Apply[f,association,{Depth[association]-2}]

End[];


(* ::Input::Initialization:: *)
HessianRoot::usage="HessianRoot[S,t,\[Tau]] calculates the Hessian root \!\(\*SqrtBox[FractionBox[SuperscriptBox[\((2  \[Pi])\), \(2\)], \(\*SuperscriptBox[\(\[ImaginaryI]\), \(2\)] Det[\*SubsuperscriptBox[\(\[PartialD]\), \({t, tt}\), \(2\)]S]\)]]\).";
Begin["`Private`"];
HessianRoot[S_,t_,\[Tau]_]:=Sqrt[(2\[Pi])/(I Derivative[0,2][S][t,t-\[Tau]])] Sqrt[(2\[Pi] Derivative[0,2][S][t,t-\[Tau]])/(I( Derivative[2,0][S][t,t-\[Tau]]Derivative[0,2][S][t,t-\[Tau]]-Derivative[1,1][S][t,t-\[Tau]]^2))]
End[];


(* ::Input::Initialization:: *)
FindStokesTransitions::usage="FindStokesTransitions[S,\[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]\[LeftAssociation]1\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(11\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(11\)]\)},2\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(12\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(12\)]\)}\[RightAssociation],\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(2\)]\)\[RightArrow]\[LeftAssociation]1\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(21\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(21\)]\)},2\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(22\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(22\)]\)}\[RightAssociation],\[Ellipsis]\[RightAssociation]] finds the set {{\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(S\)]\)},{\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(AS\)]\)},n} of the Stokes and anti-Stokes transition energies for the given set of saddle points, where Re(S) changes sign after the \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(S\)]\) and Im(S) changes sign after the \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(AS\)]\), and n is the index of the member of the pair that should be chosen after the transition (taken as the member with a positive imaginary part of the action at the largest \!\(\*SubscriptBox[\(\[CapitalOmega]\), \(i\)]\) in the given keys).

FindStokesTransitions[S,\[LeftAssociation]\!\(\*SubscriptBox[\(label\), \(1\)]\)\[RightArrow]\[LeftAssociation]\!\(\*SubscriptBox[\(\[CapitalOmega]\), \(1\)]\)\[RightArrow]\[Ellipsis]\[RightAssociation]\[RightAssociation]] finds the Stokes transitions for the given set of saddle-point curve pairs, and returns them labeled with the \!\(\*SubscriptBox[\(label\), \(i\)]\).";
FindStokesTransitions::saddleno="FindStokesTransitions called with `1` of `2` saddle-point sets of length different from 2, with set length structure `3`. Excluding those sets from the calculation.";
FindStokesTransitions::multipleS="FindStokesTransitions found multiple Stokes transitions; using `1` to return a single transition.";
FindStokesTransitions::multipleAS="FindStokesTransitions found multiple anti-Stokes transitions; using `1` to return a single transition.";
ChooserFunction::usage="ChooserFunction is an option for FindStokesTransitions that specifies which transition to take if there are multiple transitions in the given dataset. The default is Last and gives the one with higher energy; to get the full set of transitions found use Full or Identity.";
ReperiodingFunction::usage="ReperiodingFunction is an option for FindStokesTransitions, SPAdipole and UAdipole which specifies a function f[t,\[Tau]] of recombination time t and excursion time \[Tau] that will be used to re-period the pairs {t,\[Tau]} into the form {t+f[t,\[Tau]],\[Tau]}. The default is Function[0], but if pairs are split it can be useful to set ReperiodingFunction to Function[{t,\[Tau]},Floor[-Re[t-\[Tau]],\!\(\*FractionBox[\(2  \[Pi]\), \(\[Omega]\)]\)]] for \[Omega] the carrier frequency. In general, however, it is preferable to do this in a single go using ReperiodSaddles.";

Begin["`Private`"];

Protect[ReperiodingFunction,ChooserFunction];
Options[FindStokesTransitions]={ReperiodingFunction->Function[{t,\[Tau]},0],ChooserFunction->Automatic};

FindStokesTransitions[S_,deeperAssociation_/;Depth[deeperAssociation]==5,options:OptionsPattern[]]:=Map[
FindStokesTransitions[S,#,options]&,
deeperAssociation
]

FindStokesTransitions[S_,saddlesAssociation_,options:OptionsPattern[]]:=Block[{reducedSaddlesAssociation,actionList,signsList,s,processor},
reducedSaddlesAssociation=KeySort[Select[saddlesAssociation,Length[#]==2&]];
If[Length[saddlesAssociation]-Length[reducedSaddlesAssociation]>0,
Message[FindStokesTransitions::saddleno,Length[saddlesAssociation]-Length[reducedSaddlesAssociation],Length[saddlesAssociation],First/@Tally/@Split[Values[Length/@saddlesAssociation]]
]
];
actionList=ReIm[
Map[(*reduces each \[CapitalOmega]\[Rule]\[LeftAssociation]1\[Rule]Subscript[S, 1],2\[Rule]Subscript[S, 2]\[RightAssociation] to \[CapitalOmega]\[RightArrow](Subscript[S, 1]-Subscript[S, 2])*)
Apply[Subtract],
MapIndexed[(*reduces each \[CapitalOmega]\[Rule]\[LeftAssociation]1\[Rule]{Subscript[t, 1],Subscript[\[Tau], 1]},2\[Rule]{Subscript[t, 2],Subscript[\[Tau], 2]}\[RightAssociation] to \[CapitalOmega]\[Rule]\[LeftAssociation]1\[Rule]Subscript[S, 1],2\[Rule]Subscript[S, 2]\[RightAssociation]*)
With[{t=#1[[1]]+OptionValue[ReperiodingFunction][#1[[1]],#2[[2]]],\[Tau]=#1[[2]],\[CapitalOmega]=#2[[1,1]]},
S[t,t-\[Tau]]-\[CapitalOmega] t
]&
,reducedSaddlesAssociation,{2}
]]
];
signsList=Sign[Times[
Rest[actionList],
AssociationThread[Rest[Keys[actionList]],Most[Values[actionList]]]
]];
processor=OptionValue[ChooserFunction]/.{Automatic->Last,Full->Identity};
If[Length[Keys[Select[signsList,#[[1]]<0&]]]>1,Message[FindStokesTransitions::multipleS,processor]];
If[Length[Keys[Select[signsList,#[[2]]<0&]]]>1,Message[FindStokesTransitions::multipleAS,processor]];
{
processor[Keys[Select[signsList,#[[1]]<0&]]/.{{}->{Missing["No transition"]}}],
processor[Keys[Select[signsList,#[[2]]<0&]]/.{{}->{Missing["No transition"]}}],
Sign[Last[actionList][[2]]]/.{1->2,-1->1}
}
]

End[];


(* ::Input::Initialization:: *)
SPAdipole::usage="SPAdipole[S,prefactor,\[CapitalOmega],{t,\[Tau]}] returns the saddle-point approximation amplitude corresponding to action S[t,t-\[Tau]]-\[CapitalOmega]t and the given prefactor[t,t-\[Tau]].

SPAdipole[S,prefactor,\[CapitalOmega],\[LeftAssociation]1\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)},2\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(2\)]\)},\[Ellipsis]\[RightAssociation]] returns the total harmonic-dipole contribution in the saddle-point approximation from the specified saddle points.

SPAdipole[S,prefactor,\[CapitalOmega],\[LeftAssociation]1\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)},2\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(2\)]\)}\[RightAssociation],transition] uses the given Stokes transition set to drop the relevant saddle after the anti-Stokes transition.";
SPAdipole::wrongno="SPAdipole called with a Stokes transition but with an input association of length `1` at harmonic energy \[CapitalOmega]=`2`. Reverting to unstructured evaluation.";
SPAdipole::invldtrns="SPAdipole called with invalid Stokes transition set `1`. Reverting to unstructured evaluation.";

Begin["`Private`"];

Options[SPAdipole]={ReperiodingFunction->Function[{t,\[Tau]},0]};

SPAdipole[S_,prefactor_,\[CapitalOmega]_,{t_,\[Tau]_},options:OptionsPattern[]]:=Block[{tr=t+OptionValue[ReperiodingFunction][t,\[Tau]]},
HessianRoot[S,tr,\[Tau]]prefactor[tr,tr-\[Tau]]Exp[-I S[tr,tr-\[Tau]]+I \[CapitalOmega] tr]
]
SPAdipole[S_,prefactor_,\[CapitalOmega]_,times_Association,options:OptionsPattern[]]:=Block[{},
Total[SPAdipole[S,prefactor,\[CapitalOmega],#,options]&/@times]
]

SPAdipole[S_,prefactor_,\[CapitalOmega]_,times_Association,transition_,options:OptionsPattern[]]:=Block[{},
If[!NumberQ[transition[[2]]],Message[SPAdipole::invldtrns,transition];Return[SPAdipole[S,prefactor,\[CapitalOmega],times] ]];
If[Length[times]!=2,Message[SPAdipole::wrongno,Length[times],\[CapitalOmega]];Return[SPAdipole[S,prefactor,\[CapitalOmega],times,options]]];
If[\[CapitalOmega]<transition[[1]],
SPAdipole[S,prefactor,\[CapitalOmega],times,options],
SPAdipole[S,prefactor,\[CapitalOmega],KeySelect[times,#==transition[[3]]&],options]
]
]

End[];


(* ::Input::Initialization:: *)
UAdipole::usage="UAdipole[S,prefactor,\[CapitalOmega],\[LeftAssociation]1\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(1\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(1\)]\)},2\[RightArrow]{\!\(\*SubscriptBox[\(t\), \(2\)]\),\!\(\*SubscriptBox[\(\[Tau]\), \(2\)]\)},\[Ellipsis]\[RightAssociation],transition] returns the total harmonic-dipole contribution in the uniform approximation from the specified saddle points, using the action S[t,t-\[Tau]]-\[CapitalOmega]t and prefactor[t,t-\[Tau]], and taking the given Stokes transition set as a reference.";
UAdipole::saddleno="UAdipole called with `1` time pairs at \[CapitalOmega]=`2`. Reverting to the saddle-point approximation for this set.";
UAdipole::invldtrns="UAdipole called with invalid Stokes transition set `1`. Reverting to the saddle-point approximation for this set.";

Begin["`Private`"];

Options[UAdipole]={ReperiodingFunction->Function[{t,\[Tau]},0]};

UAdipole[S_,prefactor_,\[CapitalOmega]_,times_,transition_,options:OptionsPattern[]]:=(
If[Length[times]!= 2,Message[UAdipole::saddleno,Length[times],\[CapitalOmega]];Return[SPAdipole[S,prefactor,\[CapitalOmega],times] ]];
If[!NumberQ[transition[[2]]],Message[UAdipole::invldtrns,transition];Return[SPAdipole[S,prefactor,\[CapitalOmega],times] ]];
Block[
{A1,A2,S1,S2,Ss,Sm,z,
t1=times[1][[1]]+OptionValue[ReperiodingFunction][times[1][[1]],times[1][[2]]],\[Tau]1=times[1][[2]],
t2=times[2][[1]]+OptionValue[ReperiodingFunction][times[2][[1]],times[2][[2]]],\[Tau]2=times[2][[2]]},
A1=HessianRoot[S,t1,\[Tau]1]prefactor[t1,t1-\[Tau]1];
S1=S[t1,t1-\[Tau]1]- \[CapitalOmega] t1;
A2=HessianRoot[S,t2,\[Tau]2]prefactor[t2,t2-\[Tau]2];
S2=S[t2,t2-\[Tau]2]- \[CapitalOmega] t2;
Ss=(S1+S2)/2;Sm=(S1-S2)/2;
If[\[CapitalOmega]<transition[[2]],z=(-(3/2)Sm)^(2/3),z=(-(3/2)Sm)^(2/3) Exp[I (transition[[3]]/.{2->-1,1->1}) (2 \[Pi])/3]];
Sqrt[6\[Pi] Sm]Exp[-I Ss+I \[Pi]/4]((A1-I A2)/2 AiryAi[-z]/Sqrt[z]+I (A1+I A2)/2 AiryAi'[-z]/z)
])

End[];


(* ::Input::Initialization:: *)
HCAdipole::usage="HCAdipole[cutoff,prefactor,k,\[CapitalOmega]] calculates the HCA dipole from a cutoff (in the format \[LeftAssociation]\"t\"\[Rule]\[Placeholder],\"\[Tau]\"\[Rule]\[Placeholder],\"S\"\[Rule]\[Placeholder],\"\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)S\"\[Rule]\[Placeholder],\"\!\(\*SubsuperscriptBox[\(d\), \(t\), \(3\)]\)S\"\[Rule]\[Placeholder],\"\!\(\*SubsuperscriptBox[\(\[PartialD]\), \(tt\), \(2\)]\)S\"\[Rule]\[Placeholder]\[RightAssociation]) and prefactor (as a function prefactor[t,tt]) using the branch-cut factor \!\(\*SuperscriptBox[\(\[ExponentialE]\), \(\[ImaginaryI]\\\ 2  \[Pi]\\\ k/3\)]\) at harmonic frequency \[CapitalOmega].";
HCAdipole::wrongcutoff="The cutoff `1` does not contain the right set of keys.";

Begin["`Private`"];
HCAdipole[cutoff_,prefactor_,k_,\[CapitalOmega]_,options:OptionsPattern[]]:=Block[{thc,tthc,\[Tau]hc,S0,\[CapitalOmega]hc,d2Sdtt2,A},
If[Sort[Keys[cutoff]]!=Sort[{"t","\[Tau]","S","\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)S","\!\(\*SubsuperscriptBox[\(d\), \(t\), \(3\)]\)S","\!\(\*SubsuperscriptBox[\(\[PartialD]\), \(tt\), \(2\)]\)S"}],Message[HCAdipole::wrongcutoff,cutoff];Return[]];
thc=cutoff["t"];
tthc=cutoff["t"]-cutoff["\[Tau]"];
\[Tau]hc=cutoff["\[Tau]"];
S0=cutoff["S"];
\[CapitalOmega]hc=cutoff["\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)S"];
A=cutoff["\!\(\*SubsuperscriptBox[\(d\), \(t\), \(3\)]\)S"]/2;
d2Sdtt2=cutoff["\!\(\*SubsuperscriptBox[\(\[PartialD]\), \(tt\), \(2\)]\)S"];
prefactor[thc,tthc]Sqrt[(2\[Pi])/(I d2Sdtt2)] (2\[Pi])/(E^(I 2\[Pi] k/3) A^(1/3)) Exp[-I S0+I \[CapitalOmega] thc]AiryAi[(\[CapitalOmega]hc-\[CapitalOmega])/(E^(I 2\[Pi] k/3) A^(1/3))]
]
End[];


(* ::Input::Initialization:: *)
EndPackage[];


(* ::Input::Initialization:: *)
DistributeDefinitions["RBSFA`"];