Graphene.mpl

#!/home/tzdyrski/bin/maple -q
#
# Solves the Van der Pol oscillator using the method of multiple scales.
#
# Written for Maple 2018
# Thomas Zdyrski (tzdyrski@ucsd.edu)
# 14 April 2018

## Define PDE
#PDE :=
#  C*D[1,1,1](n1)(r,t)
#  +A*D[2](n1)(r,t)
#  +B*n1(r,t)*D[1](n1)(r,t)
#  +Disp*D[1,1](n1)(r,t):
#
### Remove dissipation
##eval(%, Disp=0):
#
### Use nice values
##parameters := {A=1,B=-1,C=1,Disp=-1}:
#
## Use realistic values
#parameters := {A=2.0,B=-0.83,C=1.4,Disp=-0.27}:
#
### Add steady, dissipationless ICs matching derivation
##ICs := {n1(r,0) = signum(B)*sech(sqrt(abs(B)/12/C)*r)^2,
##  n1(-10,t)=n1(10,t), D[1](n1)(-10,t)=D[1](n1)(10,t), D[1,1](n1)(-10,t)=D[1,1](n1)(10,t)}:
#
## Add steady, dissipationless ICs normalized to unit width
#ICs := {n1(r,0) = 12*C/B*sech(r)^2,
#  n1(-10,t)=n1(10,t), D[1](n1)(-10,t)=D[1](n1)(10,t), D[1,1](n1)(-10,t)=D[1,1](n1)(10,t)}:
#
## Eval PDE and ICs with parameters
#PDE := eval(PDE, parameters):
#ICs := eval(ICs, parameters):
#
## Solve PDE
#pds:=pdsolve(PDE, ICs, numeric, spacestep=1/256, timestep=1/60):
#
## Output to file
#plotsetup(ps,plotoutput="plot.ps",plotoptions="color"):
#
## Plot 3d mesh
#pds:-plot3d(t=0..1,r=-3..3,axes=boxed):
#plots[display](%);
#
## Generate and display plots at different times
#p1:=pds:-plot(t=0.0, numpoints=100, color=red):
#plots[display](p1);
#p2:=pds:-plot(t=0.2, numpoints=100, color=blue):
#plots[display](p2);
#p3:=pds:-plot(t=0.4, numpoints=100, color=green):
#plots[display](p3);
#p4:=pds:-plot(t=0.6, numpoints=100, color=purple):
#plots[display](p4);
#
#
## Display different times on a single plot
#plots[display]({p1,p2,p3,p4}, title=`Density profile at t=0,0.2,0.4,0.6`);
#
#
## Define number of points per unit length
#ppu:=5:
#
## Generate function n and sample entropy divergence ($\partial_{\mu}s^{\mu} = (\partial_x n)^2$)
#N := subs(pds:-value(output = listprocedure), n1(r,t)):
#entropy_div := ((r,t) -> (256*(N(r+1/125,t) - N(r-1/125,t)))^2)/350:
#int_entropy_div := (t -> add(entropy_div(1/(ppu+0.01)*j, t),
#  j=ppu*(-10)..ppu*(10))):
#
## Plot entropy divergence vs x and t surface plot
#A := Array([seq([seq(entropy_div(1/(ppu+0.01)*j, (1/20)*k), j =
#  ppu*(-10)..ppu*10)], k = 0 ..  20)],datatype=float[8]):
#plots[surfdata](A, axes=box);
#
## Plot density and entropy divergence at different times
#plot_list := []:
#for time_t in [0.1,0.6,1.2] do
#  plot(entropy_div(r,time_t),r=-10..10, color="green"):
#  plot_list := [op(plot_list), %]:
#  plot(N(r,time_t),r=-10..10, color="red"):
#  plot_list := [op(plot_list), %]:
#end do:
#
#plots[display](plot_list);
#
## Plot integrated entropy divergence
#plot(int_entropy_div(t),t=0.1..2):
#plots[display](%);
#
#quit:

#with(Physics):
#
#Setup(mathematicalnotation=true,
#  dimension=4,
#  metric=`+++-`,
#  spaceindices=lowercaselatin,quiet):
#
## Define our system of equations
#Coordinates(X=cartesian,quiet):
#
## Define Gamma
##Gamma := unapply(1/sqrt(1-u(X)^2),X):
#
## Define other functions
#{U[~mu] = Gamma(X)*Vector([u(X),0,0,1]),
#A[~mu]=[0,0,0,phi(X)]}:
#Define(op(%)):
#
## AA:=4*Pi*ee^2*d1*d2/kappa/(d1+d2):
## hbar = v_F = k_e = k_B = 1
#
## Define F after A is defined (since it involves derivatives)
#F[mu,nu]=d_[mu](A[nu]) - d_[nu](A[mu]):
#Define(%):
#
#Proj[mu,nu]=g_[mu,nu]+U[mu].U[nu]:
#Define(%):
#
##{J[~mu]=-ee*n(X)*U[~mu],
##T[~mu,~nu]=(e(X) + P(X))*U[~mu]*U[~nu] + P(X)*g_[~mu,~nu]}:
#{J[~mu]=-ee*n(X)*U[~mu] + 1/ee*sigma_Q*Proj[~mu,~nu].(temp(X)*d_[nu]((chempot_on_temp(X)))
#  + ee*F[nu,rho].U[~rho]),
#T[~mu,~nu]=(e(X) + P(X))*U[~mu]*U[~nu] + P(X)*g_[~mu,~nu] -
#  eta.Proj[~mu,~rho].Proj[~nu,~sigma].(d_[rho](U[sigma]) +
#  d_[sigma](U[rho]) -2/d*g_[rho,sigma].d_[alpha](U[~alpha])) -
#  zeta*Proj[~mu,~nu].d_[alpha](U[~alpha])}:
#Define(op(%)):
#
#Define( S[~mu]=[s(X),0,0,0]):
#
## Define equations
#PDEs := [
#d_[mu](J[~mu])/(-ee)=0,
#d_[mu](T[~0,~mu]) = F[~0,~mu].J[mu],
#d_[mu](T[~i,~mu]) = F[~i,~mu].J[mu],
#d_[mu](S[~mu](X))=-(J[~nu]+ee*n(X)*U[~nu]).(-1/ee*temp(X)*d_[nu](chempot_on_temp(X))-F[nu,alpha].U[~alpha])/temp(X) - (T[~mu,~nu]-(e(X) + P(X))*U[~mu]*U[~nu] - P(X)*g_[~mu,~nu]).d_[nu](U[mu])/temp(X)
#]:
#
## Sum over indices
#PDEs := SumOverRepeatedIndices(%):
#
## Re-write momentum equation using energy equation
#eval(PDEs[3], ~i=1):
#Expand(% - u(X)*PDEs[2]):
#
## Re-combine into PDE system
#PDEs := [op([1..2],PDEs),%,PDEs[4]]:
#
## Nondimensionalize (LHS dimensional, RHS dimensionless*dimensional_quantity)
#PDEs := eval(convert(PDEs,diff), {
#  n(X) = N__dim*n(X),
#  u(X) = U__dim*u(X),
#  e(X) = E__dim*e(X),
#  P(X) = P__dim*P(X),
#  chempot_on_temp(X) = Mu__dim/T__dim*chempot_on_temp(X),
#  temp(X) = T__dim*temp(X),
#  AA = AA__dim*AA, # O(AA) = O(d_i)*O(A=4\pi/\kappa from notes); this notation matches write-up
#  phi(X) = Phi__dim*phi(X),
#  ':-diff' = ((f,dt) -> Diff__dim*':-diff'(f,dt)),
#  s(X) = S__dim*s(X),
#  sigma_Q = Sigma_Q__dim*sigma_Q,
#  eta = Eta__dim*eta,
#  zeta = Zeta__dim*zeta
#}):
#
## Remove common factors
#PDEs := expand([
#  PDEs[1]/(Diff__dim*N__dim*U__dim), # divide by order(partial_x) order(n) order(u)
#  PDEs[2]/(Diff__dim*P__dim*U__dim), # divide by order(partial_x) order(P) order(u)
#  PDEs[3]/(Diff__dim*P__dim*U__dim), # divide by order(partial_x) order(P) order(u)
#  PDEs[4]/(Diff__dim*S__dim/epsilon) # divide by order(s)*order(partial_x)/epsilon ; Note: we want an extra factor of epsilon multiplying partial_x (s) since order(partial_x s) = epsilon order(partial_x) order(s), but all the other terms are quadratic, so for eg. order( (partial_x n)^2 ) = epsilon^2 order(partial_x)^2 order(n)^2
#  ]):
#
### Define new parameters for convience
## epsilon = delta n / n (aka perturbation order)
## epsilon^m = (Mu__dim/T__dim)^2
##  => epsilon^(abs(m)) = max((Mu__dim/T__dim)^2,(T__dim,Mu__dim)^2)
## epsilon^p = Diff__dim*Eta__dim/P__dim/epsilon
## epsilon^q = T__dim*epsilon^p*epsilon^(abs(m))*Sigma_Q__dim/Eta__dim
#
## Make choices (with k__B=hbar=v__F=l__ref=1)
#PDEs := eval['recurse'](%, {
#  Mu__dim = epsilon^(m/2)*T__dim, # definition of epsilon^m
#  Phi__dim = N__dim*AA__dim, # Electrostatic calculation: O(phi) = O(n)*O(AA) = O(n)*O(d_i)*O(A=4\pi/\kappa from notes)
#  E__dim = P__dim, # scale-invariant requirement
#  P__dim = (T__dim*Mu__dim)^((d+1)/2)*epsilon^(-abs(m)*((d+1)/4)), # Thermo requirement #1
#  N__dim = T__dim^((d-1)/2)*Mu__dim^((d+1)/2)*epsilon^(-abs(m)*((d-1)/4)), # Thermo requirement #2
#  D__dim = epsilon^(1/2)/Diff__dim, # diff(n,x) and d1*d2*diff(n,x$3) must differ by epsilon
#  AA__dim = P__dim/N__dim^2, # nondipersive electromagnetic term must enter at leading order
#  U__dim = 1, # 1st order u0=0 equations
#  Diff__dim = P__dim/Eta__dim*epsilon^(p+1), # Definition of epsilon^p
#  T__dim = epsilon^(q/2-p/2+abs(m)/2)*sqrt(Eta__dim/Sigma_Q__dim), # Definition of epsilon^q
#  Sigma_Q__dim = epsilon^(sigma_Q_order), # Definition of sigma_Q_order
#  Eta__dim = epsilon^(eta_order), # Definition of eta_order
#  Zeta__dim = epsilon^(zeta_order), # Definition of zeta_order
#  S__dim = epsilon^(eta_order+1+q*d/2+p*(-d+2)/2+m*(d+1)/4+abs(m)*(d-1)/4+eta_order*(d-2)/2-sigma_Q_order*d/2-1+2-p+min(p,p+zeta_order-eta_order,q,q+m/2+abs(m)/2,q+m+abs(m))) # order(s) = order(eta*diff/T*epsilon^(-1+2-p+min(p,p+zeta_order-eta_order,q,q+m/2+abs(m)/2,q+m+abs(m)); see writeup, and note that the -1 comes from the fact that order(partial_x s) = epsilon*order(partial)*order(s), with order(partial_x s) given in the writeup
#}):
#
#alias('X'=X):
#
## Make u, n, P, e, Gamma, chempot, and temp independent of x2, x3
#Gamma := unapply(GGamma(x1,x4),X):
#u := unapply(uu(x1,x4),X):
#n := unapply(nn(x1,x4),X):
#P := unapply(PP(x1,x4),X):
#e := unapply(EE(x1,x4),X):
#phi := unapply(Phi(x1,x4),X):
#chempot_on_temp := unapply(Chempot_on_temp(x1,x4),X):
#temp := unapply(Temp(x1,x4),X):
#s := unapply(ss(x1,x4),X):
#PDEs := PDEs:
#unassign('Gamma','u','n','P','e','phi','chempot_on_temp','temp','s'):
#
## Replace temporary names
#PDEs :=
#subs({GGamma=Gamma,uu=u,nn=n,PP=P,EE=e,Phi=phi,Chempot_on_temp=chempot_on_temp,Temp=temp,ss=s},%):
#
## Simplify using assumptions
#PDEs := simplify(PDEs) assuming d>0,epsilon>0,m::real,q>=0,p>=0,eta_order>=0,zeta_order>=0,sigma_Q_order>=0:
#print(PDEs):
#lprint(PDEs):
#quit:

Graphene := module()

export
  UnitTests,
  Graphene:

UnitTests := proc(
{
quiet_tests::truefalse:=false, # display test results
## Specify factors of epsilon (see writeup for meaning of m, p, q, etc)
nondim_orders::set:={m=1,p=0,q=0,sigma_Q_order=1/2,eta_order=0,zeta_order=0},
## Adiabatic or isothermal conditions
adiabatic::truefalse:=false,
# Frame choice
v0equals0::truefalse:=true, # choose u0 such that v0:=0
v1equals0::truefalse:=false, # choose U1 such that v1:=0
v1HalfEquals0::truefalse:=false, # choose U1 such that it is half the required value for v1:=0 (this can give a nicer plot)
## Printing options
print_num_values::truefalse:=false, # Print numerical values of variables
print_dimensional::truefalse:=false, # Print dimensional form of variables
print_coeffs_of_dim_vars::truefalse:=false, # Print KdVB coeffs using dimensional variables
print_observables_without_free_param::truefalse:=false, # Print observables after eliminating free param c1
verify_full_sol::truefalse:=false # Verify final form of fullly substituted solution
}
)
local
#  PhysicsProcTest,
#  MetricBugTest,
  pdes,
  delta__def,
  indicator,
  deltas,
  clean_pdes,
  entropy_div,
  calculated_entropy_div,
  min_order,
  adiabatic_local,
  comb_result,
  comb_clean,
  DiracEq,
  FermiTerms,
  result_assume,
  coefficients,
  KdVB,
  compat_assume,
  coefficients_adi,
  coefficients_iso,
  sigma_q_terms,
  result_gen,
  KdVB__deriv,
  U2,
  HOT,
  compat_gen,
  result_HOT,
  cleaner_HOT,
  sol,
  physical_units,
  LL,
  kappa__norm,
  nondim_vals,
  thermo_coeffs,
  dim_to_nondim,
  physical_constants,
  nondim_to_dim,
  char_units,
  common_factor:

description "Run unit tests.":

## Test bug fixes
#PhysicsProcTest := proc()
#  # Should return x1
#  # Without bug fix, returns X[1]
#  uses Physics:
#
#  Coordinates(X,quiet):
#
#  return X[1]:
#end proc:
#
#CodeTools:-Test(
#  PhysicsProcTest(),
#  x1,
#  label="Physics inside proc Bug Fix Test"):
#
#MetricBugTest :=  proc()
#  # Should return -diff(n(X),x1)*u(X),
#  # Without bug fix, returns -diff(n(X),x1)*u(X)-2*diff(n(X),x4)
#  uses Physics:
#
#  Setup(metric=`+++-`,quiet):
#
#  # Define our system of equations
#  Coordinates(Y, quiet):
#
#  # Define tensors
#  {U[~mu] = [u(X),0,0,1],
#  A[~mu]=[0,0,0,n(X)]}:
#  Define(op(%), quiet):
#
#  F[mu,nu]=d_[mu](A[nu]) - d_[nu](A[mu]):
#  Define(%, quiet):
#
#  # Multiply tensors
#  TensorArray(F[~0,~mu].U[mu]):
#
#  return %:
#end proc:
#
#CodeTools:-Test(
#  MetricBugTest(),
#  -diff(n(X),x1)*u(X),
#  label="Physics metric Bug Fix Test"):

Graphene(2, 't_num'=2, ':-nondim_orders'=nondim_orders):
pdes := %:

# Define a Kronecker delta function
# Source: https://stackoverflow.com/questions/15389034/is-there-a-kronecker-symbol-in-maple
delta__def := delta = ((m,n)-> `if`(evalb(m < n)::truefalse, # check if the two inputs can be compared
  `if`(m=n,1,0), # if they can be compared, check for equality
  'procname'(m,n))): # otherwise, return delta(m,n)

# Define indicator function mapping true->1 and false->0
indicator := x -> eval(x, [true=1,false=0]):

# Define deltas as a list so they can be applied at will, while retaining the option to *not* apply them
deltas :=
eval['recurse']({
  P1_corrections_delta = delta(m,-1), # If in the Fermi-regime with m=-1, the thermodynamics requires P__1 depends on (T0/mu0)^2, which changes additional equations
  Thermo_sigma_Q_delta = Heaviside(-m), # If in the Fermi-regime, the sigma_Q*chempot and sigma_Q*temp terms are the same order as the sigma_Q*AA term and must be included
  T1_corrections_delta = delta(m,-1), # If in the adiabatic Fermi-regime with m=-1, the thermodynamics requires T__1 depends on P__2 and n__2, which changes the way equations are eliminated
  sigma_Q_delta = eval(sigma_Q*delta(q,0), nondim_orders), # Define nondim_orders dependent dissipative parameters
  eta_delta = eval(eta*delta(p,0),nondim_orders),
  zeta_delta = eval(zeta*delta(p,0),nondim_orders)
}, nondim_orders union {delta__def}) union {delta__def}:

# Remove extra factors of Gamma
pdes := map(x -> [x[1]/Gamma,x[2],x[3]*Gamma,x[4]], pdes):
pdes := eval(pdes, Gamma = 1/sqrt(1-u__0(x,t__0,t__1)^2)):

# Equilibrium Solution
pdes := eval(pdes, {
  n__0 = ((x,t0,t1) -> n0),
  P__0 = ((x,t0,t1) -> P0),
  e__0 = ((x,t0,t1) -> e0),
  u__0 = ((x,t0,t1) -> u0),
  chempot__0 = ((x,t0,t1) -> mu0),
  temp__0 = ((x,t0,t1) -> T0),
  s__0 = ((x,t0,t1) -> s0)
}):

# Check form of equations
clean_pdes :=
[[
diff(n__1(x,t__0,t__1),t__0)+Gamma^2*n0*u0*diff(u__1(x,t__0,t__1),t__0)+u0*diff(n__1(x,t__0,t__1),x)+n0*Gamma^2*diff(u__1(x,t__0,t__1),x)=0,
Gamma^2*diff(e__1(x,t__0,t__1),t__0)+Gamma^2*u0^2*diff(P__1(x,t__0,t__1),t__0)+2*u0*(e0+P0)*Gamma^4*diff(u__1(x,t__0,t__1),t__0)+(1+u0^2)*(e0+P0)*Gamma^4*diff(u__1(x,t__0,t__1),x)+u0*Gamma^2*diff(e__1(x,t__0,t__1)+P__1(x,t__0,t__1),x)+AA*n0*u0*Gamma^2*diff(n__1(x,t__0,t__1),x)+AA*n0^2*u0^2*Gamma^4*diff(u__1(x,t__0,t__1),x)=0,
Gamma^3*(e0+P0)*diff(u__1(x,t__0,t__1),t__0)+Gamma*u0*diff(P__1(x,t__0,t__1),t__0)+u0*Gamma^3*(e0+P0)*diff(u__1(x,t__0,t__1),x)+Gamma*diff(P__1(x,t__0,t__1),x)+AA*n0*Gamma*diff(n__1(x,t__0,t__1),x)+AA*n0^2*u0*Gamma^3*diff(u__1(x,t__0,t__1),x)=0
],
[
diff(n__2(x,t__0,t__1),t__0)+Gamma^2*n0*u0*diff(u__2(x,t__0,t__1),t__0)+u0*diff(n__2(x,t__0,t__1),x)+n0*Gamma^2*diff(u__2(x,t__0,t__1),x)
=-diff(n__1(x,t__0,t__1),t__1)-n0*u0*Gamma^2*diff(u__1(x,t__0,t__1),t__1)-u0*u__1(x,t__0,t__1)*Gamma^2*diff(n__1(x,t__0,t__1),t__0)-Gamma^2*(Gamma^2*(1+2*u0^2)*n0*u__1(x,t__0,t__1)+u0*n__1(x,t__0,t__1))*diff(u__1(x,t__0,t__1),t__0)-Gamma^2*(u0*Gamma^2*(2+u0^2)*n0*u__1(x,t__0,t__1)+n__1(x,t__0,t__1)+n0*u0*u__1(x,t__0,t__1))*diff(u__1(x,t__0,t__1),x)-Gamma^2*u__1(x,t__0,t__1)*diff(n__1(x,t__0,t__1),x)+Gamma*u0*AA*sigma_Q_delta/ee^2*diff(n__1(x,t__0,t__1),t__0,x)+Gamma*AA*sigma_Q_delta/ee^2*diff(n__1(x,t__0,t__1),x,x)+Gamma^3*u0*AA*sigma_Q_delta/ee^2*n0*(diff(u__1(x,t__0,t__1),x,x)+u0*diff(u__1(x,t__0,t__1),t__0,x))+Thermo_sigma_Q_delta*Gamma*sigma_Q_delta/ee^2*(u0^2*diff(chempot__1(x,t__0,t__1),t__0,t__0)+2*u0*diff(chempot__1(x,t__0,t__1),t__0,x)+diff(chempot__1(x,t__0,t__1),x,x))-Thermo_sigma_Q_delta*Gamma*sigma_Q_delta/ee^2*mu0/T0*(u0^2*diff(temp__1(x,t__0,t__1),t__0,t__0)+2*u0*diff(temp__1(x,t__0,t__1),t__0,x)+diff(temp__1(x,t__0,t__1),x,x)),
Gamma^2*diff(e__2(x,t__0,t__1),t__0)+Gamma^2*u0^2*diff(P__2(x,t__0,t__1),t__0)+2*u0*(e0+P0)*Gamma^4*diff(u__2(x,t__0,t__1),t__0)+(1+u0^2)*(e0+P0)*Gamma^4*diff(u__2(x,t__0,t__1),x)+u0*Gamma^2*diff(e__2(x,t__0,t__1)+P__2(x,t__0,t__1),x)+AA*n0*u0*Gamma^2*diff(n__2(x,t__0,t__1),x)+AA*n0^2*u0^2*Gamma^4*diff(u__2(x,t__0,t__1),x)
=-2*(e0+P0)*u0*Gamma^4*diff(u__1(x,t__0,t__1),t__1)-Gamma^2*diff(e__1(x,t__0,t__1),t__1)-Gamma^2*u0^2*diff(P__1(x,t__0,t__1),t__1)-2*(e0+P0)*Gamma^6*(1+3*u0^2)*u__1(x,t__0,t__1)*diff(u__1(x,t__0,t__1),t__0)-2*(e__1(x,t__0,t__1)+P__1(x,t__0,t__1))*Gamma^4*u0*diff(u__1(x,t__0,t__1),t__0)-2*u0*u__1(x,t__0,t__1)*Gamma^4*diff(e__1(x,t__0,t__1)+P__1(x,t__0,t__1),t__0)-2*(3+u0^2)*u0*u__1(x,t__0,t__1)*Gamma^6*(e0+P0)*diff(u__1(x,t__0,t__1),x)-Gamma^4*(1+u0^2)*(e__1(x,t__0,t__1)+P__1(x,t__0,t__1))*diff(u__1(x,t__0,t__1),x)-(1+u0^2)*u__1(x,t__0,t__1)*Gamma^4*diff(e__1(x,t__0,t__1)+P__1(x,t__0,t__1),x)+Gamma^5*u0*(zeta_delta+2*eta_delta*(1-1/d))*(u0^2*diff(u__1(x,t__0,t__1),t__0,t__0)+2*u0*diff(u__1(x,t__0,t__1),t__0,x)+diff(u__1(x,t__0,t__1),x,x))-AA*n0*u__1(x,t__0,t__1)*Gamma^4*diff(n__1(x,t__0,t__1),x)-AA*u0*n__1(x,t__0,t__1)*Gamma^2*diff(n__1(x,t__0,t__1),x)-AA*n0*u0*Gamma^2*d1*d2/3*diff(n__1(x,t__0,t__1),x,x,x)-2*AA*n0^2*u0*(1+u0^2)*u__1(x,t__0,t__1)*Gamma^6*diff(u__1(x,t__0,t__1),x)-2*AA*n0*u0^2*n__1(x,t__0,t__1)*Gamma^4*diff(u__1(x,t__0,t__1),x)-AA*n0^2*u0^2*Gamma^4*d1*d2/3*diff(u__1(x,t__0,t__1),x,x,x)-AA*n0*u0^2*u__1(x,t__0,t__1)*Gamma^4*diff(n__1(x,t__0,t__1),x),
Gamma^3*(e0+P0)*diff(u__2(x,t__0,t__1),t__0)+Gamma*u0*diff(P__2(x,t__0,t__1),t__0)+u0*Gamma^3*(e0+P0)*diff(u__2(x,t__0,t__1),x)+Gamma*diff(P__2(x,t__0,t__1),x)+AA*n0*Gamma*diff(n__2(x,t__0,t__1),x)+AA*n0^2*u0*Gamma^3*diff(u__2(x,t__0,t__1),x)
=-(e0+P0)*Gamma^3*diff(u__1(x,t__0,t__1),t__1)-u0*Gamma*diff(P__1(x,t__0,t__1),t__1)-Gamma^3*(2*u0*u__1(x,t__0,t__1)*Gamma^2*(e0+P0)+(e__1(x,t__0,t__1)+P__1(x,t__0,t__1)))*diff(u__1(x,t__0,t__1),t__0)-u__1(x,t__0,t__1)*Gamma*diff(P__1(x,t__0,t__1),t__0)-Gamma^3*(u0*(e__1(x,t__0,t__1)+P__1(x,t__0,t__1))+(1+u0^2)*u__1(x,t__0,t__1)*Gamma^2*(e0+P0))*diff(u__1(x,t__0,t__1),x)-AA*n__1(x,t__0,t__1)*Gamma*diff(n__1(x,t__0,t__1),x)-AA*n0*Gamma*d1*d2/3*diff(n__1(x,t__0,t__1),x,x,x)-AA*n0^2*(1+u0^2)*u__1(x,t__0,t__1)*Gamma^5*diff(u__1(x,t__0,t__1),x)-AA*n0^2*u0*Gamma^3*d1*d2/3*diff(u__1(x,t__0,t__1),x,x,x)-2*AA*n0*u0*n__1(x,t__0,t__1)*Gamma^3*diff(u__1(x,t__0,t__1),x)+Gamma^4*(zeta_delta+2*eta_delta*(1-1/d))*(u0^2*diff(u__1(x,t__0,t__1),t__0,t__0)+2*u0*diff(u__1(x,t__0,t__1),t__0,x)+diff(u__1(x,t__0,t__1),x,x))
]]:

clean_pdes := eval(clean_pdes, ee=1):

[pdes[1,1..3],pdes[2,1..3]] - clean_pdes:

eval['recurse'](%, deltas): # Apply deltas
convert(%,diff):
eval(%, Gamma=1/sqrt(1-u0^2)):
map2(map, x->lhs(x)-rhs(x), %):
radnormal(%):
if quiet_tests then
        CodeTools:-Test(%, [[0,0,0],[0,0,0]], label="Check Form of Equations", quiet=true):
else
        CodeTools:-Test(%, [[0,0,0],[0,0,0]], label="Check Form of Equations"):
end if:

# Separate entropy divergence from pdes
entropy_div := [pdes[1,4],pdes[2,4]]:
if adiabatic then
  pdes := [pdes[1,1..3],pdes[2,1..3]]:
else
  pdes := [pdes[1,[1,3]],pdes[2,[1,3]]]:
end if:

## Check entropy divergence
eval(entropy_div[2], Gamma = 1/sqrt(1-u0^2)):
calculated_entropy_div := eval(%,u0=0):

# Define useful combination of parameters (see nondim_orders of entropy divergence)
min_order := min(p,p+zeta_order-eta_order,q,q+m/2+abs(m)/2,q+m+abs(m)):

# Known form of entropy divergence for u0=0
D[1](s__1)(x,t__0,t__1) =
((D[1](chempot__1)(x,t__0,t__1)-mu0/T0*D[1](temp__1)(x,t__0,t__1))*delta(q+m+abs(m),min_order)
  +AA*D[1](n__1)(x,t__0,t__1)*delta(q,min_order))^2*sigma_Q/T0
+1/T0*D[1](u__1)(x,t__0,t__1)^2*(zeta*delta(p+zeta_order-eta_order,min_order)+2*eta*(1-1/d)*delta(p,min_order)):

# Compare calculated and known forms
%-calculated_entropy_div:

eval(%, nondim_orders): # Apply nondim_orders
eval(%, delta__def): # Apply definition of delta

verify(lhs(%)-rhs(%),0,simplify):
if quiet_tests then
        CodeTools:-Test(%, true, label="Check Entropy Divergence", quiet=true):
else
        CodeTools:-Test(%, true, label="Check Entropy Divergence"):
end if:

## Check cleaner form of first-order solutions

## Define thermodynamic expansion of pressure P in terms of mu and T
## For Dirac (m>0):
## P = T_0^(d+1)*(C__0+eps^m*C__1*(mu/T)^2+eps^(2*m)*C__2*(mu/T)^4)
## For Fermi (m<0):
## P = abs(mu_0)^(d+1)*(C__0+eps^(-m)*C__1*(T/mu)^2+eps^(-2*m)*C__2*(T/mu)^4)
## For Fermi, m=-1, we have
#C__0 := 8*surfNsphere__dMinus1/(2*Pi*hbar*v__F)*(1/(2*(d+1)*d)):
#C__1 := 8*surfNsphere__dMinus1/(2*Pi*hbar*v__F)*(Pi^2/12):
#C__2 := 8*surfNsphere__dMinus1/(2*Pi*hbar*v__F)*(7*Pi^4/720):

# Thermodynamic Relations
if adiabatic then
  eval['recurse'](pdes, {
    e0 = P0*d,
    e__1 =       ((x,t0,t1) -> P__1(x,t0,t1)*d),
    chempot__1 = ((x,t0,t1) -> mu0*(n__1(x,t0,t1)/n0/d*(1+delta(m,-1)*(d-1)) - P__1(x,t0,t1)/P0*(d-1)/(d+1)*delta(m,-1))*Heaviside(-m)
                              +mu0*(n__1(x,t0,t1)/n0 - P__1(x,t0,t1)/P0*(d-1)/(d+1) + (mu0/T0)^2*(C__1/C__0*(d-1)/(d+1)-2*C__2/C__1)*delta(m,1))*Heaviside(m)),
    temp__1 =    ((x,t0,t1) -> T0*(-(d-1)*chempot__1(x,t0,t1)/mu0+d/2*C__0/C__1*(mu0/T0)^2*(P__2(x,t0,t1)/P0-n__2(x,t0,t1)/n0*(d+1)/d-chempot__1(x,t0,t1)^2/mu0^2*(d+1)/2-3/d*C__2/C__0*(T0/mu0)^4))*delta(m,-1)*Heaviside(-m)
                              +T0/(d+1)*(P__1(x,t0,t1)/P0-C__1/C__0*(mu0/T0)^2*delta(m,1))*Heaviside(m)),
    e__2 =       ((x,t0,t1) -> P__2(x,t0,t1)*d)
  }):
else
  eval['recurse'](pdes, {
    e0 = P0*d,
    e__1 =       ((x,t0,t1) -> P__1(x,t0,t1)*d),
    e__2 =       ((x,t0,t1) -> P__2(x,t0,t1)*d),
    temp__1 =    ((x,t0,t1) -> T1),
    chempot__1 = ((x,t0,t1) -> P0/n0*(d+1)/d*(n__1(x,t0,t1)/n0 - C__1/C__0*(d-1)/(d+1)*(T0/mu0)^2*delta(m,-1))*Heaviside(-m)
                              +mu0*(n__1(x,t0,t1)/n0 - 2*C__2/C__1*(mu0/T0)^2*delta(m,1))*Heaviside(m)),
    P__1 =       ((x,t0,t1) -> P0*(n__1(x,t0,t1)/n0*(d+1)/d + C__1/C__0/d*(T0/mu0)^2*delta(m,-1))*Heaviside(-m)
                              +P0*(C__1/C__0*(mu0/T0)^2*delta(m,1))*Heaviside(m)),
    P__2 =       ((x,t0,t1) -> P0*(n__2(x,t0,t1)/n0*(d+1)/d + n__1(x,t0,t1)^2/n0^2*(d+1)/2/d^2 + C__1/C__0*(d-1)/d^2*(T0/mu0)^2*(n__1(x,t0,t1)/n0 - 3/2*C__1/C__0*(d-1)/(d+1)*(T0/mu0)^2)*delta(m,-1) + 3*C__2/C__0/d*(T0/mu0)^4*delta(m,-1) + 1/d*C__1/C__0*(T0/mu0)^2*delta(m,-2))*Heaviside(-m)
                              +P0*(2*n__1(x,t0,t1)/n0*C__1/C__0*(mu0/T0)^2*delta(m,1) - 3*C__2/C__0*(mu0/T0)^4*delta(m,1) + C__1/C__0*(mu0/T0)^2*delta(m,2))*Heaviside(m))
  }):
  if eval(m,nondim_orders)>0 then
    eval(%, mu0 = n0/2/C__1/T0^(d-1)):
  end if:
end if:

eval(%, nondim_orders): # Apply nondim_orders

# Move to lhs
pdes := map2(map,x -> lhs(x) - rhs(x) = 0, %):

# For some reason, maple evaluates arguments like adiabatic in the set
# definitions, but doesn't evaluate local variables.
# Also, it doesn't like the eval['recurse'](pdes[1]),%) statement when
# the argument adiabatic is already evaluated to true/false, so replace
# the argument with a local variable to delay evluation
adiabatic_local := adiabatic:

{Gamma = 1/sqrt(1-u0^2),
K0 = (d+1)/d*indicator(adiabatic_local or eval(m,nondim_orders)<0),
v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic_local)
        + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic_local),
n__1 = ((x,t__0,t__1) -> n1(x+v0*t__0,t__1)+_F3(t__1)+_F4(x-u0*t__0,t__1)*indicator(adiabatic_local)),
P__1 = ((x,t__0,t__1) -> P0*(d+1)/d/n0*n1(x+v0*t__0,t__1)+_F2(t__1)-AA*n0*Gamma^2*_F4(x-u0*t__0,t__1)), # only applied for adiabatic_local since, for isothermal, P__1 substitution was already made
u__1 = ((x,t__0,t__1) -> -1/n0/Gamma^2*(v0+u0)/(1+u0*v0)*n1(x+v0*t__0,t__1)+_F1(t__1))}:

eval['recurse'](pdes[1],%):
eval(%, nondim_orders): # Apply nondim_orders
eval(%, Gamma = 1/sqrt(1-u0^2)):
expand(%):
simplify(%):

convert(%,set):
if quiet_tests then
        CodeTools:-Test(%,{0=0},label="Check cleaner form of first-order solutions", quiet=true):
else
        CodeTools:-Test(%,{0=0},label="Check cleaner form of first-order solutions"):
end if:

## Check combination of first-order equations
pdes[1]:

if adiabatic then
  Gamma^2*AA*n0*d*diff(%[1],x) +diff(%[2],x)+u0*diff(%[2],t__0)
    -Gamma*((u0^2+d)*diff(%[3],t__0) + u0*(d+1)*diff(%[3],x)):
  diff(%, t__0) + u0*diff(%,x):
else
  Gamma^2*d*(AA*n0*diff(%[1],x)+P0/n0*K0*(u0*diff(%[1],t__0)+diff(%[1],x)))-Gamma*d*(diff(%[2],t__0)+u0*diff(%[2],x)):
end if:
lhs(%)-rhs(%):
comb_result := %:

# Check cleaner form
u__1(x,t__0,t__1): # place this here for convience
if adiabatic then
  Gamma^2*(d+1)*P0*(u0^2-d)*diff(%, t__0$2)-2*Gamma^2*(d+1)*P0*u0*(d-1)*diff(%,x,t__0)+(AA*n0^2*d+Gamma^2*(d+1)*P0*(1-d*u0^2))*diff(%, x$2):
  Gamma^2*(diff(%, t__0) + u0*diff(%, x)):
else
 Gamma^2*d*((AA*n0^2+Gamma^2*P0*(K0-u0^2*(d+1)))*diff(%,x$2)+2*Gamma^2*P0*u0*(K0-(d+1))*diff(%,x,t__0)+Gamma^2*P0*(K0*u0^2-(d+1))*diff(%,t__0$2)):
end if:
comb_clean := %:

comb_result - comb_clean:

eval(%, Gamma=1/sqrt(1-u0^2)):
eval(%, K0 = (d+1)/d*Heaviside(-m)):
eval(%, nondim_orders):
radnormal(%):
expand(%):
simplify(%):

if quiet_tests then
        CodeTools:-Test(%,0,label="Check combination of first-order equations", quiet=true):
else
        CodeTools:-Test(%,0,label="Check combination of first-order equations"):
end if:

## Verify cleaner form is equivalent (assuming n__2,u__2,P__2 right moving waves at speed v0)

# Only consider third-order equation
pdes[2]:

# Combine equations assuming travelling waves
if adiabatic then
  DiracEq:=
  Gamma^2*AA*n0*d*%[1]+%[2]*(1+u0*v0)
    -Gamma*%[3]*(v0*(u0^2+d) + u0*(d+1)):

  pdes[2]:
  FermiTerms:=-%[2]*(1+u0*v0)+Gamma*%[3]*(v0*(u0^2+d)+u0*(d+1))
    +AA*d*n0^2/P0/(d+1)*(Gamma*u0*%[3]-%[2]):

  (u0+v0)*DiracEq+T1_corrections_delta*Gamma/2*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*(1+u0*v0)^2*diff(FermiTerms,x):
else
  Gamma^2*d*(AA*n0+P0/n0*K0*(1+u0*v0))*%[1]-Gamma*d*(u0+v0)*%[2]:
  DiracEq:=eval(%, K0 = (d+1)/d*Heaviside(-m)):

  pdes[2]:
  FermiTerms:=0:

  DiracEq:
end if:

eval['recurse'](%, deltas): # Apply deltas

# Travelling wave solutions
eval(%, {
  n__2 = ((x,t0,t1) -> n2(x+v0*t0,t1)),
  u__2 = ((x,t0,t1) -> u2(x+v0*t0,t1)),
  P__2 = ((x,t0,t1) -> P2(x+v0*t0,t1)) # only applied for adiabatic since, for isothermal, P__2 substitution was already made
}):

# Insert first-order solutions
eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic),
  n__1 = ((x,t__0,t__1) -> n1(x+v0*t__0,t__1)),
  P__1 = ((x,t__0,t__1) -> P0*(d+1)/d/n0*n1(x+v0*t__0,t__1)+P1_corrections_delta*P0*(C__1/C__0/d*(T0/mu0)^2)), # only applied for adiabatic since, for isothermal, P__1 substitution was already made
  u__1 = ((x,t__0,t__1) -> -1/n0/Gamma^2*(v0+u0)/(1+u0*v0)*n1(x+v0*t__0,t__1)+U1)}
):
eval['recurse'](%, deltas): # Apply deltas
result_assume := eval(%, nondim_orders): # Apply nondim_orders

# Replace argument with left-moving coordinate
v0*t__0+x=X:
eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):
eval(%, nondim_orders): # Apply nondim_orders
result_assume := subs(%, result_assume):

## Cleaner form of KdV
if adiabatic then
  coefficients:=
  {
  A = 2*Gamma^2*(u0+v0)/(1+u0*v0)*P0*(d+1)/n0*(v0*(d-u0^2)+u0*(d-1)),
  B = -Gamma^2*P0*(d+1)/n0^2*(u0+v0)/(1+u0*v0)*(d-1)/d*(3*d*(u0+v0)^2-(1+u0*v0)^2),
  C = -AA*n0*d*d1*d2/3*(u0+v0)/(1+u0*v0),
  F = Gamma^2*P0*(d+1)/n0*(u0+v0)/(1+u0*v0)*(2*U1*Gamma^2*(d-1)*(u0+v0)*(1+u0*v0) +P1_corrections_delta/d*C__1/C__0*(T0/mu0)^2*(d*(u0+v0)^2-(1+u0*v0)^2)),
  G = -d*Gamma*(1+u0*v0)/n0*((n0^2*AA^2/(1+u0*v0)+mu0*AA*n0/d*Gamma^2*(1+u0*v0)*Heaviside(-m-1.5))*sigma_Q_delta+Gamma^2*(u0+v0)^2*(zeta_delta+2*eta_delta*(1-1/d)))
  }:
else
  coefficients:=
  {
  A = 2*Gamma^2*(u0+v0)/(1+u0*v0)*P0*d/n0*(v0*(d+1-u0^2*K0)+u0*(d+1-K0)),
  B = -Gamma^2*P0/n0^2/d*(u0+v0)/(1+u0*v0)*(d^2*(u0+v0)^2*(4*(d+1)-K0*(d+3))+(1+u0*v0)^2*((d+1)*Heaviside(-m)-K0*d^2)), # Heaviside comes from P__2 = Heaviside*n1^2 term
  C = -AA*n0*d*d1*d2/3*(u0+v0)/(1+u0*v0),
  F = Gamma^2*d*P0/n0*(u0+v0)/(1+u0*v0)*(2*Gamma^2*(u0+v0)*(1+u0*v0)*(d+1-K0)*U1
      +C__1/C__0*(mu0/T0)^(2*m)*(   (d+1)*(u0+v0)^2*(1/d*delta(m,-1)+delta(m,1))
                                                   -(1+u0*v0)^2*((d-1)/d^2*delta(m,-1)+2*delta(m,1)))),
  G = -Gamma^3*(1+u0*v0)/n0*(Gamma^2*(u0+v0)*(1+u0*v0)*(d+1)*sigma_Q_delta*P0^2/n0^2*((v0*(d+1-K0*u0^2)+u0*(d+1-K0))*(d*(u0+v0)^2/(1+u0*v0)^2+Thermo_sigma_Q_delta-K0*d/(d+1)))
    +d*(u0+v0)^2*(zeta_delta+2*eta_delta*(1-1/d)))
  }:
end if:

# Define KdV-Burgers eq
KdVB := -(A*D[2](n1)(X,t__1)+F*D[1](n1)(X,t__1)+B*n1(X,t__1)*D[1](n1)(X,t__1) +C*D[1,1,1](n1)(X,t__1) + G*D[1,1](n1)(X,t__1)):
# Don't add higher order terms to KdVB to get full (n1,u1,etc as well as n2,u2,etc
# terms), cleaner compat equations: these terms drop out automatically
# since we assumed u2,n2,etc were right moving waves at the correct
# phase speed
if adiabatic then
  0=(u0+v0)*KdVB-(1+u0*v0)^2*T1_corrections_delta/2*Gamma*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*diff(Eval(KdVB,sigma_Q=0),X):
else
  0=KdVB:
end if:

eval(%, eval(coefficients, nondim_orders)):
eval['recurse'](%, deltas): # Apply deltas
if eval(m,nondim_orders)>0 then # Replace mu0 using thermo
  eval(%, mu0 = n0/2/C__1/T0^(d-1)):
end if:
convert(%,D):
value(%): # Now that we've applied the form of G, evaluate sigma_Q=0

compat_assume := %:

# Compare to cleaner equations
result_assume - compat_assume:
lhs(%)-rhs(%)=0:

eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):
eval(%, nondim_orders): # Apply nondim_orders
eval(%, Gamma = 1/sqrt(1-u0^2)):

radnormal(%):
expand(%):
simplify(%):
if quiet_tests then
        CodeTools:-Test(%, 0=0, label="Check cleaner form of KdV (assuming right-moving)", quiet=true):
else
        CodeTools:-Test(%, 0=0, label="Check cleaner form of KdV (assuming right-moving)"):
end if:

## Check relation between adiabatic and isothermal coefficients

coefficients_adi :=
[
  A = 2*Gamma^2*(u0+v0)/(1+u0*v0)*P0*(d+1)/n0*(v0*(d-u0^2)+u0*(d-1)),
  B = -Gamma^2*P0*(d+1)/n0^2*(u0+v0)/(1+u0*v0)*(d-1)/d*(3*d*(u0+v0)^2-(1+u0*v0)^2),
  C = -AA*n0*d*d1*d2/3*(u0+v0)/(1+u0*v0),
  F = Gamma^2*P0*(d+1)/n0*(u0+v0)/(1+u0*v0)*(2*U1*Gamma^2*(d-1)*(u0+v0)*(1+u0*v0) +P1_corrections_delta/d*C__1/C__0*(T0/mu0)^2*(d*(u0+v0)^2-Correction*(1+u0*v0)^2)),
  G = -d*Gamma*(1+u0*v0)/n0*(n0^2*AA^2/(1+u0*v0)*sigma_Q_delta+Gamma^2*(u0+v0)^2*(zeta_delta+2*eta_delta*(1-1/d)))
]:
eval(%,{sigma_Q_delta=sigma_Q,eta_delta=eta,zeta_delta=zeta,P1_corrections_delta=delta(m,-1)}): # Apply fake deltas
coefficients_adi := map(x->lhs(x)-rhs(x),%):
coefficients_adi := eval(%, v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))):

coefficients_iso :=
[
  A = 2*Gamma^2*(u0+v0)/(1+u0*v0)*P0*d/n0*(v0*(d+1-u0^2*K0)+u0*(d+1-K0)),
  B = -Gamma^2*P0/n0^2/d*(u0+v0)/(1+u0*v0)*(d^2*(u0+v0)^2*(4*(d+1)-K0*(d+3))+(1+u0*v0)^2*((d+1)*Heaviside(-m)-K0*d^2)), # Heaviside comes from P__2 = Heaviside*n1^2 term
  C = -AA*n0*d*d1*d2/3*(u0+v0)/(1+u0*v0),
  F = Gamma^2*d*P0/n0*(u0+v0)/(1+u0*v0)*(2*Gamma^2*(u0+v0)*(1+u0*v0)*(d+1-K0)*U1
      +C__1/C__0*(mu0/T0)^(2*m)*(   (d+1)*(u0+v0)^2*(1/d*delta(m,-1)+delta(m,1))
                                                   -(1+u0*v0)^2*((d-1)/d^2*delta(m,-1)+2*delta(m,1)))),
  G = -Gamma^3*(1+u0*v0)/n0*(Gamma^2*(u0+v0)*(1+u0*v0)*(d+1)*sigma_Q_delta*P0^2/n0^2*((v0*(d+1-K0*u0^2)+u0*(d+1-K0))*(d*(u0+v0)^2/(1+u0*v0)^2+Thermo_sigma_Q_delta-K0*d/(d+1)))
    +d*(u0+v0)^2*(zeta_delta+2*eta_delta*(1-1/d)))
]:
eval(%,{sigma_Q_delta=sigma_Q,eta_delta=eta,zeta_delta=zeta,Thermo_sigma_Q_delta=delta(m,-1)}): # Apply fake deltas
coefficients_iso := map(x->lhs(x)-rhs(x),%):
coefficients_iso := eval(%, v0 = (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)):
coefficients_iso := eval(%, K0 = (d+1)/d*Heaviside(-m)):

# Apply Dirac m=1 to G's sigma_Q terms
eval(coefficients_iso[5], sigma_Q=0):
coefficients_iso[5]-%:
sigma_q_terms := eval(%, m=1):

# Apply Fermi m=-1 to G's eta/zeta (ie, all other) terms
eval(coefficients_iso[5], sigma_Q=0):
eval(%, m=-1):
% + sigma_q_terms:
coefficients_iso := [op(1..4,coefficients_iso), %]:

coefficients_iso := eval(coefficients_iso, m=-1): # Apply Fermi m=-1 to all other coefficients

coefficients_adi := eval(coefficients_adi, m=-1): # Apply Fermi m=-1 (\mathcal{F} term doesn't match for m=1 Dirac)

# Compare
coefficients_adi - coefficients_iso:

eval(%, delta__def): # Apply definition of delta
eval(%, Correction = (d-1)/(d+1)):
eval(%, Gamma=1/sqrt(1-u0^2)):
radnormal(%):
expand(%):
simplify(%) assuming d>1, u0^2<1:
convert(%, set):

CodeTools:-Test(%, {0}, label="Check relation between adiabatic and isothermal coefficients", quiet=quiet_tests):

# Verify cleaner form is equivalent (w/o assuming n__2,u__2,P__2 right moving waves at speed v0)

# Only consider third-order equation
pdes[2]:

# Combine equations assuming travelling waves
if adiabatic then
  DiracEq:=
  Gamma^2*AA*n0*d*diff(%[1],x) +diff(%[2],x)+u0*diff(%[2],t__0)
    -Gamma*((u0^2+d)*diff(%[3],t__0) + u0*(d+1)*diff(%[3],x)):

  pdes[2]:
  FermiTerms:=-(diff(%[2],x)+u0*diff(%[2],t__0))+Gamma*((u0^2+d)*diff(%[3],t__0)+u0*(d+1)*diff(%[3],x))
    +AA*d*n0^2/P0/(d+1)*diff(Gamma*u0*%[3]-%[2],x):

  diff(DiracEq,t__0) + u0*diff(DiracEq,x)
    +T1_corrections_delta*Gamma/2*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*(diff(FermiTerms,x$2)
    +2*u0*diff(FermiTerms,x,t__0)+u0^2*diff(FermiTerms,t__0$2)):
else
  Gamma^2*d*(AA*n0*diff(%[1],x)+P0/n0*K0*(u0*diff(%[1],t__0)+diff(%[1],x)))-Gamma*d*(diff(%[2],t__0)+u0*diff(%[2],x)):
  DiracEq:=eval(%, K0 = (d+1)/d*Heaviside(-m)):

  pdes[2]:
  FermiTerms:=0:

  DiracEq:
end if:
eval['recurse'](%, deltas): # Apply deltas

# Insert first-order solutions
eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic),
  n__1 = ((x,t__0,t__1) -> n1(x+v0*t__0,t__1)),
  P__1 = ((x,t__0,t__1) -> P0*(d+1)/d/n0*n1(x+v0*t__0,t__1)+P1_corrections_delta*P0*(C__1/C__0/d*(T0/mu0)^2)), # only applied for adiabatic since, for isothermal, P__1 substitution was already made
  u__1 = ((x,t__0,t__1) -> -1/n0/Gamma^2*(v0+u0)/(1+u0*v0)*n1(x+v0*t__0,t__1)+U1)}
):
eval['recurse'](%, deltas): # Apply deltas
result_gen := eval(%, nondim_orders): # Apply nondim_orders

# Replace argument with left-moving coordinate
v0*t__0+x=X:
eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):
eval(%, nondim_orders): # Apply nondim_orders
result_gen := subs(%, result_gen):

## Cleaner form of KdV
# Define KdV-Burgers eq
KdVB := -(A*D[2](n1)(X,t__1)+F*D[1](n1)(X,t__1)+B*n1(X,t__1)*D[1](n1)(X,t__1) +C*D[1,1,1](n1)(X,t__1) + G*D[1,1](n1)(X,t__1)):
KdVB__deriv := diff(KdVB,X):
# Define for conventience
U2 := u__2(x,t__0,t__1):
# Define eqn for higher order terms
HOT := -Gamma^4*P0*d*(d+1-u0^2*K0)*(v0p*v0m*diff(U2,x$2)-(v0p+v0m)*diff(U2,x,t__0)+diff(U2,t__0$2)):
# Add higher order terms to KdVB to get full (n1,u1,etc as well as n2,u2,etc
# terms), cleaner compat equations
if adiabatic then
  diff(HOT,t__0)+u0*diff(HOT,x)-T1_corrections_delta/2*Gamma*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*(diff(HOT,x$2)+2*u0*diff(HOT,x,t__0)+u0^2*diff(HOT,t__0$2))
    =(u0+v0)*diff(KdVB__deriv,X)-(1+u0*v0)^2*T1_corrections_delta/2*Gamma*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*Eval(diff(KdVB__deriv,X$2),sigma_Q=0):
else
  HOT=KdVB__deriv:
end if:

eval(%, eval(coefficients, nondim_orders)):
eval['recurse'](%, deltas): # Apply deltas
if eval(m,nondim_orders)>0 then # Replace mu0 using thermo
  eval(%, mu0 = n0/2/C__1/T0^(d-1)):
end if:
convert(%,D):
value(%): # Now that we've applied the form of G, evaluate sigma_Q=0

compat_gen := %:

# Compare to cleaner equations
result_gen - compat_gen:
lhs(%)-rhs(%)=0:

# We need the definitions for both left and right moving waves
eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0p = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic),
  v0m = (-u0*(d-1)/(d-u0^2)-sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))-1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):

eval['recurse'](%, {
  K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):
eval(%, nondim_orders): # Apply nondim_orders
eval(%, Gamma = 1/sqrt(1-u0^2)):

radnormal(%):
expand(%):

simplify(%):

if quiet_tests then
        CodeTools:-Test(%, 0=0, label="Check cleaner form of KdV (w/o assuming right-moving)", quiet=true):
else
        CodeTools:-Test(%, 0=0, label="Check cleaner form of KdV (w/o assuming right-moving)"):
end if:

## Check diff eq for HOT in terms of left and right moving coordinates, chi_0^(+) and chi_0^(-)
if adiabatic then
  diff(HOT,t__0)+u0*diff(HOT,x)-T1_corrections_delta/2*Gamma*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*(diff(HOT,x$2)+2*u0*diff(HOT,x,t__0)+u0^2*diff(HOT,t__0$2)):
else
  HOT:
end if:

eval['recurse'](%, deltas): # Apply deltas
result_HOT := %:

if adiabatic then
  Gamma^4*P0*d*(d+1-u0^2*K0)*(v0p-v0m)^2*diff(u__2(chi0p,chi0m,t__1),chi0p,chi0m):
  (u0+v0p)*diff(%,chi0p)+(u0+v0m)*diff(%,chi0m)
    -T1_corrections_delta/2*Gamma*C__0/C__1*sigma_Q_delta*mu0^3/T0^2*(d+1)/n0*((1+u0*v0p)^2*diff(%,chi0p,chi0p)+2*(1+u0*v0p)*(1+u0*v0m)*diff(%,chi0p,chi0m)+(1+u0*v0m)^2*diff(%,chi0m,chi0m)):
else
  Gamma^4*P0*d*(d+1-u0^2*K0)*(v0p-v0m)^2*diff(u__2(chi0p,chi0m,t__1),chi0p,chi0m):
end if:

eval['recurse'](%, deltas): # Apply deltas
PDEtools:-dchange({chi0p=x+v0p*t__0,chi0m=x+v0m*t__0}, % , [x,t__0], params={v0m,v0p}):
cleaner_HOT := %:

# Compare to cleaner equations
result_HOT - cleaner_HOT:

# We need the definitions for both left and right moving waves
eval['recurse'](%, {
  K0 = (d+1)/d*Heaviside(-m),
  v0p = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic),
  v0m = (-u0*(d-1)/(d-u0^2)-sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))-1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):

eval['recurse'](%, {
  K0 = (d+1)/d*Heaviside(-m),
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}):
eval(%, nondim_orders): # Apply nondim_orders
eval(%, Gamma = 1/sqrt(1-u0^2)):

radnormal(%):
expand(%):
if quiet_tests then
        CodeTools:-Test(%, 0, label="Check diff eq for HOT in terms of left and right moving coordinates", quiet=true):
else
        CodeTools:-Test(%, 0, label="Check diff eq for HOT in terms of left and right moving coordinates"):
end if:

## Generate solution

## General form of solution to KdV in co-moving frame
#KdV := 0 = -v*A*D(f)(x) + B*f(x)*D(f)(x) + C*D[1,1,1](f)(x):
#f(x) = 3*v*A/B*sech(sqrt(v*A/C)/2*x)^2:
#verify(pdetest(%, KdV),0);

## General form of solution to KdV in stationary frame
#KdV := 0 = A*D[2](f)(x,t)+F*D[1](f)(x,t) + B*f(x,t)*D[1](f)(x,t) + C*D[1,1,1](f)(x,t):
#f(x,t) = v*sech(sqrt(v*B/12/C)*(x - (v*B/3/A+F/A)*t))^2:
#f(x,t) = v*sech(sqrt(v*B/12/C)*(x))^2 - v/3:
#verify(pdetest(%, KdV),0);

## General form of solution to KdV Burgers
#KdVBurgers := C*D[1,1,1](n1)(X,t) +A*D[2](n1)(X,t) +B*n1(X,t)*D[1](n1)(X,t) +G*D[1,1](n1)(X,t):
#n1(X,t) = 3*G^2/25/B/C*(sech(-G/C/10*X-3*G^3/125/A/C^2*t)^2-2*tanh(-G/C/10*X-3*G^3/125/A/C^2*t)-2):
#verify(pdetest(%, KdVBurgers),0);

## Dissipationless
#sol := 3*v*A/B*sech(sqrt(v*A/C)/2*(X - (v+F/A)*t__1))^2:
#coefficients := eval(coefficients, {sigma_Q=0,eta=0,zeta=0}):

# Dissipative
# Plus
sol := 3*G^2/25/B/C*(sech(-G/C/10*X+(F*G/10/C/A+3*G^3/125/A/C^2)*t__1)^2-2*tanh(-G/C/10*X+(F*G/10/C/A+3*G^3/125/A/C^2)*t__1)+2):
## Minus
#sol := 3*G^2/25/B/C*(sech(-G/C/10*X+(F*G/10/C/A-3*G^3/125/A/C^2)*t__1)^2-2*tanh(-G/C/10*X+(F*G/10/C/A-3*G^3/125/A/C^2)*t__1)-2):

# Apply coefficients
sol := eval(sol, eval(coefficients, nondim_orders)):


## Specify numerical values of order-1, nondim variables (denoted by primes in paper); here, l__ref = hbar = v__F = k__B = 1
## Note: l__ee doesn't appear in any of the equations/definitions, so even
## after setting hbar=k__B=v__F=1, we have one undeteremined dimension (at this
## point, everything is expressed in, eg., mass units or geometrized units

## Display variables in units of (meters, seconds, kilograms, kelvin, ohms); ie, define a unit system
physical_units := {
  l__ref = LL*1e-7*Unit(`m`), # reference length 100 nm = 1e-7 m
  v__F = 1e6*Unit(`m`/`s`), # Fermi velocity 1e8 m / s
  k__B = 1.38064852e-23*Unit(`kg`*`m`^2/`s`^2/`K`), # Boltzmann constant 1.4e-19 kg m^2 / s^2 / K
  hbar = 1.0545718e-34*Unit(`kg`*`m`^2/`s`), # Planck's constant 1.055e-34 kg m^2 / s
  ee = 1.60217662e-19*Unit(`C`) # Elementary charge 1.602e-19 C
}:

# Make choice of l__ref
LL:=0.5:

# It is more useful to specify kappa unnormalized; however, we need its
# (epsilon_var) normalized value below. Note: kappa doesn't need to be
# order-unity since it doesn't appear on its own: instead, AA must be
# order-unity, and we're absorbing AA/di normalization to ensure AA is
# order-unity
kappa__norm := 1*epsilon_var^(1/2+q-sigma_Q_order):

## Note: when choosing params, we have two approaches:
# 1) set dimensional parmeters (eg T0=100K,d1=3um,etc)
#   -all observables are independent of l__ref
#   -nondimensional quantities (including c1) *do* vary with l__ref
#   -l__ref changes are dangerous b/c nondimensional parameters can be made not order-1
#   -we use this method for l__ref
# 2) set nondimensional parameters (eg n0=1,T0=1,etc)
#   -nondim parameters (including c1) are independent of l__ref
#   -dimenional quantities *do* vary with l__ref
#   -l__ref changes are dangerous b/c dimensional variables can be made unrealistic
#   -we use this method for epsilon (except for T0, where we use method 1 for both epsilon and l_ref)

# Define values of independent variables we can choose
nondim_vals :=
{
  d = 2,
  epsilon_var = 0.1,
  T0 = 60*Unit(`K`)*k__B/hbar/v__F*l__ref/epsilon_var^(q/2-p/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  n0 = 4*LL^2,
  d1 = 0.5/LL,
  d2 = 0.5/LL,
  c1 = 4*LL^d,
  kappa = kappa__norm # dielectric constant (normalized)
}:
map(x->lhs(x)=eval['recurse'](rhs(x), %),%): # evaluate nested definitions
eval(%, physical_units):
eval(%, nondim_orders): # Apply nondim_orders before deltas
combine(%, 'units'):

# In the manuscript, when enumerating the system's free parameters, we neglect
# l__ref and epsilon_var; that's because these are simply artifacts of working
# in nondimensional variables. Upon redimensionalizing, they go away.

if v0equals0 then
  % union {u0 = sqrt((1/d+AA*n0^2/P0/(d+1))/(1+AA*n0^2/P0/(d+1)))*indicator(adiabatic) # This gives v0=0
    +sqrt((K0/(d+1)+AA*n0^2/P0/(d+1))/(1+AA*n0^2/P0/(d+1)))*indicator(not adiabatic)}:
else
  % union {u0 = 0}:
end if:

if v1equals0 then
  % union {U1 = # This gives v1=0
1/6*(((-C__0*c1*d*mu0^2+3*C__1*T0^2*n0*P1_corrections_delta+C__0*c1*mu0^2)*v0^2+3*d*(C__0*c1*d*mu0^2-C__1*T0^2*n0*
P1_corrections_delta-C__0*c1*mu0^2))*u0^2+6*v0*(d-1)*(C__0*c1*d*mu0^2-P1_corrections_delta*C__1*T0^2*n0-1/3*c1*C__0*mu0^2)*u0+3*d*(
C__0*c1*d*mu0^2-C__1*T0^2*n0*P1_corrections_delta-C__0*c1*mu0^2)*v0^2-C__0*c1*d*mu0^2+3*P1_corrections_delta*C__1*T0^2*n0+c1*C__0*mu0^
2)/n0/d/C__0/mu0^2/Gamma^2/(d-1)/(v0+u0)/(u0*v0+1)*indicator(adiabatic)
+1/6*(-3*n0*(((v0+u0)^2*d^2+((-v0^2+1)*u0^2+v0^2-1)*d+(u0*v0+1)^2)*delta(m,-1)+d^2*((v0+u0)^2*d+(-2*v0^2+1)*u0^2-2*u0*v0+v0^2-2)*delta
(m,1))*C__1*(mu0/T0)^(2*m)-C__0*c1*((u0*v0+1)^2*(d+1)*Heaviside(m)+(v0+u0)^2*(K0-4)*d^3+((K0*v0^2+3*K0-4)*u0^2+8*v0*(K0-1)*u0+(3*K0-4)
*v0^2+K0)*d^2-(u0*v0+1)^2*d-(u0*v0+1)^2))/n0/d^2/C__0/Gamma^2/(v0+u0)/(u0*v0+1)/(d+1-K0)*indicator(not adiabatic)}:
elif v1HalfEquals0 then
  % union {U1 = # This is 1/2 the required value for v1=0
1/2*(
1/6*(((-C__0*c1*d*mu0^2+3*C__1*T0^2*n0*P1_corrections_delta+C__0*c1*mu0^2)*v0^2+3*d*(C__0*c1*d*mu0^2-C__1*T0^2*n0*
P1_corrections_delta-C__0*c1*mu0^2))*u0^2+6*v0*(d-1)*(C__0*c1*d*mu0^2-P1_corrections_delta*C__1*T0^2*n0-1/3*c1*C__0*mu0^2)*u0+3*d*(
C__0*c1*d*mu0^2-C__1*T0^2*n0*P1_corrections_delta-C__0*c1*mu0^2)*v0^2-C__0*c1*d*mu0^2+3*P1_corrections_delta*C__1*T0^2*n0+c1*C__0*mu0^
2)/n0/d/C__0/mu0^2/Gamma^2/(d-1)/(v0+u0)/(u0*v0+1)*indicator(adiabatic)
+1/6*(-3*n0*(((v0+u0)^2*d^2+((-v0^2+1)*u0^2+v0^2-1)*d+(u0*v0+1)^2)*delta(m,-1)+d^2*((v0+u0)^2*d+(-2*v0^2+1)*u0^2-2*u0*v0+v0^2-2)*delta
(m,1))*C__1*(mu0/T0)^(2*m)-C__0*c1*((u0*v0+1)^2*(d+1)*Heaviside(m)+(v0+u0)^2*(K0-4)*d^3+((K0*v0^2+3*K0-4)*u0^2+8*v0*(K0-1)*u0+(3*K0-4)
*v0^2+K0)*d^2-(u0*v0+1)^2*d-(u0*v0+1)^2))/n0/d^2/C__0/Gamma^2/(v0+u0)/(u0*v0+1)/(d+1-K0)*indicator(not adiabatic)
)}:
else
  % union {U1 = 0}:
end if:

nondim_vals := eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):

# Define values of independent variable set by nature; note: T0__cutoff is dependent on
# epsilon_var since T0__cutoff is a real, non-order-1, dimensionful quantity; to nondimensionalize, we use the same procedure as nondimensionalizting the temperature T0 (see paper)
{
  T0__cutoff = 1e4*Unit('K')*k__B/hbar/v__F*l__ref/epsilon_var^(q/2-p/2+abs(m)/2+eta_order/2-sigma_Q_order/2), # Temperature cutoff T0__cutoff = 10^4 kelvin
  alpha0 = 1 # bare coupling constant for graphene
}:
eval(%, physical_units):
combine(%, `units`):
eval(%, nondim_vals):
nondim_vals := nondim_vals union %: # Append to list of nondim_vals

# Calculate values of dependent variables; note: sigma_Q is a another
# dimensionless parameter, so we needed to choose its relation to epsilon_var
# (see paper for discussion)
{
  alpha = 4/(4/alpha0+ln(T0__cutoff/T0)), # renormalized (screened) electromagnetic coupling
  P0 = (abs(n0)/(d+1))^((d+1)/d)/abs(C__0)^(1/d)*signum(C__0)*Heaviside(-m)+T0^(d+1)*C__0*Heaviside(m),
  mu0 = (abs(n0)/abs(C__0)/(d+1))^(1/d)*signum(n0*C__0)*Heaviside(-m)+n0/2/C__1/T0^(d-1)*Heaviside(m),
  Gamma = 1/sqrt(1-u0^2),
  sigma_Q = 0.76/32/Pi*ln(T0__cutoff/T0)^2/epsilon_var^(sigma_Q_order), # https://doi.org/10.1103/PhysRevB.78.085416; to nondimensionalize, see paper
  eta = (0.45*T0^2/alpha^2*Heaviside(m)+1/64/Pi/alpha^2/ln(1/alpha)*mu0^2*n0/T0^2*Heaviside(-m))/epsilon_var^(eta_order), # Calculate what eta should be, just for reference
  zeta = 0/epsilon_var^(zeta_order) # Ignore this
}:
eval(%, nondim_vals):
nondim_vals := nondim_vals union %: # Append to list of nondim_vals

# Calculate variables of common groupings
{
  AA=alpha*d1*d2/kappa/(d1+d2)
}:
eval(%, nondim_vals):
nondim_vals := nondim_vals union %: # Append to list of nondim_vals

# Calculate speed (of soliton or background flow, depending on whether we chose v0=0 or u0=0)
{
  v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
          + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(-u0^2*K0+d+1)*indicator(not adiabatic)
}:
eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):
eval(%, nondim_vals):
nondim_vals := nondim_vals union %: # Append to list of nondim_vals

# Add calculated value of decay time t_d and width W
{
  t__decay = -45*A*C/4/G/B/c1,
  W = sqrt(12*C/B/c1)
 }:
eval(%, coefficients):
eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):
eval['recurse'](%, deltas): # Apply deltas
nondim_vals := nondim_vals union %: # Append to list of nondim_vals

# Add coefficient of (partial_x x)^2 in entropy divergence: div s = entropy__div__coeff*(diff(n1(x,t),x))^2
{
  entropy__div__coeff = AA^2*delta(q,min_order)*sigma_Q/T0+1/T0*v0^2/n0^2*(zeta*delta(p+zeta_order-eta_order,min_order)
    +2*eta*(1-1/d)*delta(p,min_order))
}:
eval(%, coefficients):
nondim_vals := nondim_vals union %: # Append to list of nondim_vals

# Define thermo coefficients
thermo_coeffs :=
{
  C__0 = 8*surfNsphere__dMinus1/(2*Pi)^d*(eta(0)*1/((d+1)*d)*Heaviside(-m)+(d-1)!*eta(d+1)*Heaviside(m)),
  C__1 = 8*surfNsphere__dMinus1/(2*Pi)^d*(eta(2)*Heaviside(-m)+(d-1)!*eta(d-1)/2*Heaviside(m)),
  C__2 = 8*surfNsphere__dMinus1/(2*Pi)^d*(eta(4)*(d-1)*(d-2)*Heaviside(-m)+(d-1)!*eta(d-3)/4!*Heaviside(m)),
  surfNsphere__dMinus1 = 2*Pi^(d/2)/GAMMA(d/2),
  eta = (z -> piecewise(z=1,ln(2),(1-2^(1-z))*Zeta(z)))
}:

## Calculate values of coefficients and quantities of interest
# Display -G to match notation in paper
if print_num_values then
{A_ = A, B_ = B, C_ = C, F_ = F, G_ = -G, u0_=u0, v0_ = v0, AA_ = AA, kappa_ = kappa, d1_= d1, n0_=n0,
T0_ = T0, P0_ = P0, P0__ = P0/T0^3*T0_^3, mu0_ = mu0, mu0__ = mu0/n0*n0_*T0/T0_,
sigma_Q_ = sigma_Q, sigma_Q__unnorm_ = sigma_Q*epsilon_var^(sigma_Q_order), G_/sqrt(B_*C_) = G/sqrt(B*C), alpha_ = alpha, T0__cutoff_ = T0__cutoff,
t__decay_ = t__decay, U1_=U1, c1_ = c1, eta_ = eta, eta__unnorm_ = eta*epsilon_var^(eta_order), entropy__div__coeff_ = entropy__div__coeff}:

eval(%, eval(coefficients, nondim_orders)):
eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):

eval(%, nondim_orders): # Apply nondim_orders before deltas
eval['recurse'](%, deltas): # Apply deltas
eval['recurse'](%, nondim_vals): # Give nondimensional, order-1 results
eval['recurse'](%, thermo_coeffs): # Apply thermo coeffs
eval(%, nondim_orders): # Apply nondim_orders again
eval(%, delta__def):
eval(%, d=2):
simplify(%) assuming LL>0:

#print(%):
print(evalf(%)):
end if:


## Specify values of dimesionful variables
if print_dimensional then
dim_to_nondim := {
  d__ = d,
  epsilon_var__ = epsilon_var,
  T0__dim = T0/k__B*hbar*v__F/l__ref*epsilon_var^(q/2-p/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  n0__dim = n0/l__ref^d*epsilon_var^(q*d/2-p*d/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*d/2-sigma_Q_order*d/2),
  d1__dim = d1*l__ref*epsilon_var^(-1/2-q*(d+1)/2+p*(d-1)/2-m*(d+1)/4-abs(m)*(d+1)/4-eta_order*(d-1)/2+sigma_Q_order*(d+1)/2),
  d2__dim = d2*l__ref*epsilon_var^(-1/2-q*(d+1)/2+p*(d-1)/2-m*(d+1)/4-abs(m)*(d+1)/4-eta_order*(d-1)/2+sigma_Q_order*(d+1)/2),
  kappa__dim=kappa*epsilon_var^(-1/2-q+sigma_Q_order),
  u0__dim = u0*v__F,
  U1__dim = U1*v__F*epsilon_var,
  T0__cutoff__dim = T0__cutoff/k__B*hbar*v__F/l__ref*epsilon_var^(q/2-p/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  alpha0__ = alpha0,
  alpha__ = alpha,
  P0__dim = P0*hbar*v__F/l__ref^(d+1)*epsilon_var^(q*(d+1)/2-p*(d+1)/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*(d+1)/2-sigma_Q_order*(d+1)/2),
  mu0__dim = mu0*hbar*v__F/l__ref*epsilon_var^(q/2-p/2+m/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  Gamma__ = Gamma,
  sigma_Q__dim = sigma_Q*ee^2/hbar*l__ref^(2-d)*epsilon_var^(sigma_Q_order),
  rho_Q__dim = 1/sigma_Q__dim,
  zeta__dim = zeta*hbar/l__ref^d*epsilon_var^(zeta_order),
  eta__dim = eta*hbar/l__ref^d*epsilon_var^(eta_order),
  AA__dim = AA*hbar*v__F*l__ref^(d-1)*epsilon_var^(q*(-d+1)/2+p*(d-1)/2-m*(d+1)/4-abs(m)*(d+1)/4-eta_order*(d-1)/2+sigma_Q_order*(d-1)/2),
  v0__dim = v0*v__F,
  voltage__dim = mu0__dim/ee, # gate voltage mu0/ee
  voltage__signal__dim = epsilon_var*mu0__dim/ee, # signal voltage epsilon*mu0/ee
  char__soliton__height = epsilon_var*c1*n0__dim/n0, # characteristic soliton height = eps*mag(n1) = eps*c1*order(n0*l_ref^d)/l_ref^d = eps*c1*n0__dim/n0
  K0__dim = -ee*n0__dim*u0__dim, # current density
  xi__dim = l__ref/epsilon_var^(1+q*(d+1)/2-p*(d-1)/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*(d-1)/2-sigma_Q_order*(d+1)/2), # characteristic length, inverse of partial_x
  W__dim = W*xi__dim, # width of soliton
  l__ref__dim = l__ref,
  t__ee__dim = t__ee__ratio/T0__dim, # electron-electron scattering time
  l__ee__dim = t__ee__dim*v__F, # electron-electron mean free path
  t__passing__dim = W__dim/v0__dim, # time to travel one body length
  bandwidth__dim = 1/t__passing__dim, # bandwidth of signal
  t__decay__dim = t__decay/epsilon_var*xi__dim/v__F,
  char__entropy0__dim = 1/l__ref^d*k__B*epsilon_var^(eta_order+1+q*d/2+p*(-d+2)/2+m*(d+1)/4+abs(m)*(d-1)/4+eta_order*(d-2)/2-sigma_Q_order*d/2-1+2-p+min(p,p+zeta_order-eta_order,q,q+m/2+abs(m)/2,q+m+abs(m))), # chracteristic entropy; order(s) = order(eta*diff/T*epsilon^(-1+2-p+min(p,p+zeta_order-eta_order,q,q+m/2+abs(m)/2,q+m+abs(m)))); see writeup, and note that the -1 comes from the fact that order(partial_x s) = epsilon*order(partial)*order(s), with order(partial_x s) given in the writeup
  char__entropy__div__dim = char__entropy0__dim/xi__dim*epsilon_var, # characteristic entropy density; order(partial_x s) = epsilon*order(partial)*order(s)
  char__entropy__div__int__dim = char__entropy__div__dim*xi__dim, # integral of characteristic entropy density
  graphene__specific__heat__mass__dim = graphene__specific__heat__ratio*T0__dim, # specific heat per unit mass; accurate T=0.3K-100K. Source: ln-ln interpolation from DOI:10.1103/PhysRevB.66.153408
  graphene__specific__heat__area__dim = graphene__specific__heat__mass__dim*12.01*Unit(`g`/`mol`)*6.3e-9*Unit(`mol`/`cm`^2), # specific heat per unit area
  radiative__power__dim = 4*stefan__boltzmann*graphene__emissivity*epsilon_var*(T0__dim)^4, # approximation emissivity*stefan_boltzmann*(T0__dim^4-(T0__dim+epsilon*T0__dim)^4) = 4*emissivity*boltzmann*epsilon*T0__dim^4
  graphene__thermal__cond__dim = graphene__thermal__cond__ratio*(T0__dim/Unit(`K`))^(0.557), # out-of-plane thermal conductance; accurate T=30K-100K. Source: ln-ln interpolation from DOI:10.1063/1.3245315
  joule__power__dim = K0__dim^2/sigma_Q__dim, # power density from Joule heating
  TJoule__dim = joule__power__dim/graphene__thermal__cond__dim, # Steady-state temperature increment from Joule heating (only if on substrate, ie kappa<>1)
  T__Joule__rate__dim = joule__power__dim/graphene__specific__heat__area__dim, # temperature change rate from radiation
  t__Joule__dim = epsilon_var*T0__dim/T__Joule__rate__dim, # characteristic radiation temperature change time
  T__radiat__rate__dim = radiative__power__dim/graphene__specific__heat__area__dim, # temperature change rate from radiation
  t__radiat__dim = epsilon_var*T0__dim/T__radiat__rate__dim, # characteristic radiation temperature change time
  k__F = sqrt(Pi*n0__dim), # Fermi wavenumber
  t__phonon__dim = 8*hbar^2*v__F*v__s__dim^2*graphene__density/acoustic__def__pot^2/(k__B*T0__dim)/k__F, # electron-phonon momentum relaxation time; source: DOI:10.1103/PhysRevB.76.205423
  xi__phonon__dim = t__phonon__dim*(v0__dim*u0__dim), # average distance needed for one electron to encounter one phonon
  n0__units = l__ref^(-d), # non-normalized units of n0
  d1__units = l__ref, # non-normalized units of d1
  d2__units = l__ref, # non-normalized units of d2
  AA__units = hbar*v__F*l__ref^(d-1), # non-normalized units of AA
  T0__units = 1/k__B*hbar*v__F/l__ref, # non-normalized units of T0
  P0__units = hbar*v__F/l__ref^(d+1), # non-normalized units of P0
  mu0__units = hbar*v__F/l__ref, # non-normalized units of mu0
  sigma_Q__units = ee^2/hbar*l__ref^(2-d), # non-normalized units of sigma_Q
  eta__units = hbar/l__ref^d # non-normalized units of eta
}:
dim_to_nondim := map(x->lhs(x)=eval['recurse'](rhs(x), dim_to_nondim),%): # evaluate nested definitions
dim_to_nondim := map(x->lhs(x)=eval['recurse'](rhs(x), nondim_vals),%): # Give nondimensional, order-1 results

## Define additional physical constants
physical_constants := {
  stefan__boltzmann = 5.670e-8*Unit(`W`/`m`^2/`K`^4),
  graphene__density = 7.6e-7*Unit(`kg`/`m`^2),
  t__ee__ratio = 0.1*Unit(`ps`)*100*Unit(`K`), # electron-electron scattering time at 100K, times 100K (for converting T to t__ee)
  acoustic__def__pot = 20*Unit(`eV`), # acoustic phonon deformation potential for graphene; sources vary from 10eV to 30eV
  v__s__dim = 2e4*Unit(`m`/`s`), # acoustic phonon sound speed in graphene
  graphene__thermal__cond__ratio = 5.13e6*Unit(`W`/`m`^2/`K`), # out-of-plane thermal conductivity conversion ratio; accurate T=30K-100K. Source: ln-ln interpolation from DOI:10.1063/1.3245315
  graphene__emissivity = 0.01,
  graphene__specific__heat__ratio = 1*Unit(`mJ`/`g`/`K`^2) # specific heat ratio; accurate T=0.3K-100K. Source: ln-ln interpolation from DOI:10.1103/PhysRevB.66.153408
}:

with(Units[Standard]):
dim_to_nondim:
eval(%, physical_constants):
eval(%, physical_units): # express answers in terms of (m, s, kg, K)
eval['recurse'](%, thermo_coeffs): # Apply thermo coeffs
eval(%, nondim_orders): # Apply nondim_orders
eval(%, delta__def):
eval(%, d=2):
combine(%, `units`):
simplify(%) assuming LL::positive:

#print(%):
print(evalf(%)):
unwith(Units[Standard]):

end if:

## Convert nondim, order-1 vars to dimensionful, non-order-1 vars
nondim_to_dim := {
  d = d,
  epsilon_var = epsilon_var,
  T0 = T0__dim*k__B/hbar/v__F*l__ref/epsilon_var^(q/2-p/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  n0 = n0__dim*l__ref^d/epsilon_var^(q*d/2-p*d/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*d/2-sigma_Q_order*d/2),
  d1 = d1__dim/l__ref/epsilon_var^(-1/2-q*(d+1)/2+p*(d-1)/2-m*(d+1)/4-abs(m)*(d+1)/4-eta_order*(d-1)/2+sigma_Q_order*(d+1)/2),
  d2 = d2__dim/l__ref/epsilon_var^(-1/2-q*(d+1)/2+p*(d-1)/2-m*(d+1)/4-abs(m)*(d+1)/4-eta_order*(d-1)/2+sigma_Q_order*(d+1)/2),
  kappa=kappa__dim/epsilon_var^(-1/2-q+sigma_Q_order),
  u0 = u0__dim/v__F,
  U1 = U1__dim/v__F/epsilon_var,
  T0__cutoff = T0__cutoff__dim*k__B/hbar/v__F*l__ref/epsilon_var^(q/2-p/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  alpha0 = alpha0,
  alpha = alpha,
  P0 = P0__dim/hbar/v__F*l__ref^(d+1)/epsilon_var^(q*(d+1)/2-p*(d+1)/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*(d+1)/2-sigma_Q_order*(d+1)/2),
  mu0 = mu0__dim/hbar/v__F*l__ref/epsilon_var^(q/2-p/2+m/2+abs(m)/2+eta_order/2-sigma_Q_order/2),
  Gamma = Gamma,
  sigma_Q_delta = sigma_Q__dim*hbar/l__ref^(2-d)/epsilon_var^(sigma_Q_order)*delta(q,0),
  zeta_delta = zeta__dim/hbar*l__ref^d/epsilon_var^(zeta_order)*delta(p,0),
  eta_delta = eta__dim/hbar*l__ref^d/epsilon_var^(eta_order)*delta(p,0),
  AA=AA__dim/hbar/v__F/l__ref^(d-1)/epsilon_var^(q*(-d+1)/2+p*(d-1)/2-m*(d+1)/4-abs(m)*(d+1)/4-eta_order*(d-1)/2+sigma_Q_order*(d-1)/2),
  v0 = v0__dim/v__F
}:

# Define orders of parameters (with factors of epsilon, without dimensions)
char_units :=
{
  xi__order = 1/epsilon_var^(1+q*(d+1)/2-p*(d-1)/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*(d-1)/2-sigma_Q_order*(d+1)/2), # characteristic length, inverse of partial_x
  n0__order = epsilon_var^(q*d/2-p*d/2+m*(d+1)/4+abs(m)*(d+1)/4+eta_order*d/2-sigma_Q_order*d/2)
}:

## Calculate coefficients in terms of dimensional, non-order-1 quantities
## Replace n0__order with dummy variable n0__order__ since its cleaner
if print_coeffs_of_dim_vars then
common_factor := epsilon_var^(p/2-q/2-eta_order/2+sigma_Q_order/2): # Factor out common epsilon term
{A__dim = A*(common_factor__/common_factor), B__dim =
  B*(n0__order__/n0__order)*(common_factor__/common_factor),
  C__dim = C*(xi__order/xi__order__)^2*(common_factor__/common_factor),
  F__dim = F*(common_factor__/common_factor),
  G__dim = -G*(xi__order/xi__order__)*(common_factor__/common_factor),
  v0__dim_/v_F = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)+
  (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}:

eval(%, coefficients):
#eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):

eval(%, nondim_to_dim):
eval(%, char_units):
simplify(%) assuming epsilon_var>0,l__ref>0,d>0,v__F>0,hbar>0,k__B>0,p>0,q>0,m::real:

print(%):
end if:

if print_observables_without_free_param then
## Calculate observables (there is still a free, order-1 parameter $c1$); use total $u_0' = u_0 + eps*U1$ and
# $P_0' = P_0 + eps*Theta*(T0/mu)^2...$ (see writeup)

# Note that the observables should all be independent of epsilon_var and c1
# (only dependent on their combination amplitude__dim=c1*epsilon_var*n0__char)
# Likewise, all observables should be independent of l__ref since l__ref is arbitrary
{width = sqrt(12*C/B/c1)*xi__order*l__ref, speed = v0*v__F - epsilon_var*v__F*(c1*B/3/A+F/A),
v0 = (-u0*(d-1)/(d-u0^2)+sqrt(d)/Gamma^2/(d-u0^2)*sqrt(1+AA*Gamma^2*n0^2/P0/(d+1)*(d-u0^2)))*indicator(adiabatic)
    + (u0*(K0-(d+1))+1/Gamma*sqrt(K0*(d+1)/Gamma^2+AA*n0^2/P0*(d+1-K0*u0^2)))/(d+1-u0^2*K0)*indicator(not adiabatic)
}:
#eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):

# As long as p or q is 0, then decay G is non-zero and decay_rate is finite;
# use decay_rate rather than decay_time since its easier for maple to simplify
# and cancel factors of epsilon
if eval(p,nondim_orders) = 0 or eval(q,nondim_orders) = 0 then
  % union {decay_rate = 1/(1/epsilon_var*45/4*A*C/G/B/c1*xi__order*l__ref/v__F)}:
end if:

# Write observables in terms of amplitude
eval(%, c1 = amplitude__dim/n0__order/l__ref^(-d)/epsilon_var):

eval(%, coefficients):
eval(%, nondim_to_dim): # convert to dimensionful, non-order-1 results
eval(%, char_units):
simplify(%) assuming AA__dim>0,P0__dim>0,d1__dim>0,d2__dim>0,epsilon_var>0,d::integer,d>0,
  m::real,p>=0,q>=0,sigma_Q_order>=0,eta_order>=0,zeta_order>=0,l__ref>0:

print(%):
end if:

# Verify solution
if verify_full_sol then
eval(KdVB,n1 = unapply(sol,X,t__1)):
eval(%, eval(coefficients, nondim_orders)):
eval(%, K0 = (d+1)/d*indicator(adiabatic or eval(m,nondim_orders)<0)):
eval(%, nondim_orders): # Apply nondim_orders
verify(%,0,simplify):
if quiet_tests then
        CodeTools:-Test(%, true, label="Verify Soliton Solution", quiet=true):
else
        CodeTools:-Test(%, true, label="Verify Soliton Solution"):
end if:

# Show soliton
print(`Soliton Solution`):
print(sol):
end if:
end proc:

Graphene := proc(
  e_order::integer:=1, # highest order in epsilon to keep (leading order is e_order 1)
  {leading_order_eqs::integer:=1, # Leading order of perturbation expansion for dependent variable
   nondim_orders::Or(list,set)(name=numeric):={m=1,p=0,q=0,sigma_Q_order=0,eta_order=0,zeta_order=0}, # List or set of the form name=numeric with a name for each of the following: {m,p,q,sigma_Q_order,eta_order,zeta_order} ; additionally, m must be a non-zero integer, and p and q must be non-negative integers
   t_num::integer:=e_order-leading_order_eqs+1, # number of time scales to implement
   MMS::truefalse:=false}, # if true (default: false), use the method of
                           # multiple scales instead of the Poincare-Lindstedt method
  $)
local
  entropy_div,
  PDEs,
  cleaner_PDES,
  ts,
  constraints,
  trial_functions,
  multiplier,
  j,
  orth_basis,
  tsZero,
  ICs,
  l,
  const,
  exclude_vars,
  ind_vars,
  assumptions,
  constantsToReplace,
  unit,
  pdes,
  IndDefs,
  DepDefs,
  new_ind_vars:

# Load packages
read("MultipleScales.mpl"):

interface(ansi=true, errorbreak=2):

# Define our system of equations
PDEs :=
[(-((temp(x,t)*Gamma(x,t)^2-temp(x,t))*diff(diff(chempot_on_temp(x,t),t),t)+temp(x,t)*(u(x,t)^2*Gamma(x,t)^2+1)*diff(diff(chempot_on_temp(x,t),x),x)+
2*u(x,t)*temp(x,t)*Gamma(x,t)^2*diff(diff(chempot_on_temp(x,t),t),x)+(2*u(x,t)*temp(x,t)*Gamma(x,t)*diff(Gamma(x,t),t)+2*u(x,t)^2*temp(x,t)*Gamma(x,t
)*diff(Gamma(x,t),x)+u(x,t)*Gamma(x,t)^2*diff(temp(x,t),t)+(u(x,t)^2*Gamma(x,t)^2+1)*diff(temp(x,t),x)+2*Gamma(x,t)^2*temp(x,t)*(u(x,t)*diff(u(x,t),x
)+1/2*diff(u(x,t),t)))*diff(chempot_on_temp(x,t),x)+2*diff(chempot_on_temp(x,t),t)*(temp(x,t)*Gamma(x,t)*diff(Gamma(x,t),t)+u(x,t)*temp(x,t)*Gamma(x,
t)*diff(Gamma(x,t),x)+(1/2*Gamma(x,t)^2-1/2)*diff(temp(x,t),t)+1/2*Gamma(x,t)^2*(u(x,t)*diff(temp(x,t),x)+temp(x,t)*diff(u(x,t),x))))*sigma_Q*epsilon
^(q+1/2*m+1/2*abs(m)+1)+(u(x,t)*Gamma(x,t)*epsilon^(q+1)*diff(diff(phi(x,t),t),x)*sigma_Q+Gamma(x,t)*epsilon^(q+1)*sigma_Q*diff(diff(phi(x,t),x),x)+
diff(phi(x,t),x)*sigma_Q*(u(x,t)*diff(Gamma(x,t),t)+Gamma(x,t)*diff(u(x,t),t)+diff(Gamma(x,t),x))*epsilon^(q+1)+(n(x,t)*diff(Gamma(x,t),t)+n(x,t)*
diff(Gamma(x,t),x)*u(x,t)+Gamma(x,t)*(u(x,t)*diff(n(x,t),x)+n(x,t)*diff(u(x,t),x)+diff(n(x,t),t)))*ee)*ee)/ee^2 = 0, (-zeta*d*((Gamma(x,t)^2-1)*diff(
diff(Gamma(x,t),t),t)+(2*u(x,t)*Gamma(x,t)^2-u(x,t))*diff(diff(Gamma(x,t),t),x)+(Gamma(x,t)^3-Gamma(x,t))*diff(diff(u(x,t),t),x)+diff(diff(Gamma(x,t)
,x),x)*u(x,t)^2*Gamma(x,t)^2+diff(diff(u(x,t),x),x)*u(x,t)*Gamma(x,t)^3+2*diff(Gamma(x,t),x)^2*u(x,t)^2*Gamma(x,t)+(4*u(x,t)*Gamma(x,t)*diff(Gamma(x,
t),t)+5*u(x,t)*Gamma(x,t)^2*diff(u(x,t),x)+diff(u(x,t),t)*(Gamma(x,t)^2-1))*diff(Gamma(x,t),x)+2*diff(Gamma(x,t),t)^2*Gamma(x,t)+diff(u(x,t),x)*(4*
Gamma(x,t)^2-1)*diff(Gamma(x,t),t)+diff(u(x,t),x)^2*Gamma(x,t)^3)*epsilon^(-eta_order+p+1+zeta_order)-2*(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^4+(-d*u(x,t)
^2+2*d-2)*Gamma(x,t)^2-d+1)*epsilon^(p+1)*eta*diff(diff(Gamma(x,t),t),t)-4*u(x,t)*(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^4+(-1/4*d*u(x,t)^2+5/4*d-2)*Gamma(
x,t)^2-1/4*d+1/2)*epsilon^(p+1)*eta*diff(diff(Gamma(x,t),t),x)-4*Gamma(x,t)*((1/2+(d-1/2)*u(x,t)^2)*Gamma(x,t)^4+(-1/4*d*u(x,t)^2+1/4*d-1)*Gamma(x,t)
^2-1/4*d+1/2)*epsilon^(p+1)*eta*diff(diff(u(x,t),t),x)-2*(u(x,t)^2*(u(x,t)-1)*(u(x,t)+1)*(d-1)*Gamma(x,t)^4+((3/2*d-2)*u(x,t)^2-1/2*d)*Gamma(x,t)^2+1
/2*d)*epsilon^(p+1)*eta*diff(diff(Gamma(x,t),x),x)-2*((1+(d-1)*u(x,t)^2)*Gamma(x,t)^2+d-2)*u(x,t)*Gamma(x,t)^3*epsilon^(p+1)*eta*diff(diff(u(x,t),x),
x)-2*d*u(x,t)*Gamma(x,t)^3*epsilon^(p+1)*eta*(Gamma(x,t)-1)*(Gamma(x,t)+1)*diff(diff(u(x,t),t),t)-8*(Gamma(x,t)*(u(x,t)^2*(u(x,t)-1)*(u(x,t)+1)*(d-1)
*Gamma(x,t)^2+(-1+3/4*d)*u(x,t)^2-1/4*d)*diff(Gamma(x,t),x)^2+(2*u(x,t)*Gamma(x,t)*(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^2-1/8*d*u(x,t)^2+5/8*d-1)*diff(
Gamma(x,t),t)+9/4*u(x,t)*Gamma(x,t)^2*(((d-1)*u(x,t)^2-2/9*d+7/9)*Gamma(x,t)^2+2/3*d-10/9)*diff(u(x,t),x)+2*(((d-3/8)*u(x,t)^2-1/8*d+1/8)*Gamma(x,t)^
4+(-3/16*d*u(x,t)^2+5/16*d-1/4)*Gamma(x,t)^2-1/16*d+1/8)*diff(u(x,t),t))*diff(Gamma(x,t),x)+Gamma(x,t)*(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^2+d-1/2*d*u(x
,t)^2-1)*diff(Gamma(x,t),t)^2+((((2*d-2)*u(x,t)^2+3/2-1/4*d)*Gamma(x,t)^4+(-2+3/8*d-3/8*d*u(x,t)^2)*Gamma(x,t)^2-1/8*d+1/4)*diff(u(x,t),x)+7/4*u(x,t)
*Gamma(x,t)^2*diff(u(x,t),t)*((d-2/7)*Gamma(x,t)^2-5/7*d))*diff(Gamma(x,t),t)+3/4*Gamma(x,t)^3*(((1/3+(d-1)*u(x,t)^2)*Gamma(x,t)^2+1/3*d-2/3)*diff(u(
x,t),x)^2+4/3*u(x,t)*diff(u(x,t),t)*((d-1/2)*Gamma(x,t)^2-1/4*d)*diff(u(x,t),x)+1/3*d*diff(u(x,t),t)^2*(Gamma(x,t)-1)*(Gamma(x,t)+1)))*eta*epsilon^(p
+1)+2*d*(u(x,t)*Gamma(x,t)*(e(x,t)+P(x,t))*diff(Gamma(x,t),x)+Gamma(x,t)*(e(x,t)+P(x,t))*diff(Gamma(x,t),t)+1/2*Gamma(x,t)^2*(e(x,t)+P(x,t))*diff(u(x
,t),x)+(1/2*Gamma(x,t)^2-1/2)*diff(P(x,t),t)+1/2*Gamma(x,t)^2*(diff(P(x,t),x)*u(x,t)+diff(e(x,t),x)*u(x,t)+diff(e(x,t),t))))/d = -diff(phi(x,t),x)*((
(u(x,t)^2*Gamma(x,t)^2+1)*diff(chempot_on_temp(x,t),x)+u(x,t)*Gamma(x,t)^2*diff(chempot_on_temp(x,t),t))*temp(x,t)*sigma_Q*epsilon^(q+1/2*m+1/2*abs(m
)+1)-Gamma(x,t)*ee*(u(x,t)*ee*n(x,t)+diff(phi(x,t),x)*sigma_Q*epsilon^(q+1)))/ee, (-zeta*d*((u(x,t)^2+1)*diff(diff(Gamma(x,t),t),x)+u(x,t)*diff(diff(
Gamma(x,t),t),t)+u(x,t)*diff(diff(Gamma(x,t),x),x)+diff(diff(u(x,t),t),x)*u(x,t)*Gamma(x,t)+Gamma(x,t)*diff(diff(u(x,t),x),x)+u(x,t)*Gamma(x,t)^3*
diff(u(x,t),x)^2+(Gamma(x,t)^3*diff(u(x,t),t)+(u(x,t)^2*Gamma(x,t)^2+2)*diff(Gamma(x,t),x)+u(x,t)*diff(Gamma(x,t),t)*(Gamma(x,t)^2+1))*diff(u(x,t),x)
+(u(x,t)*(Gamma(x,t)^2+1)*diff(Gamma(x,t),x)+Gamma(x,t)^2*diff(Gamma(x,t),t))*diff(u(x,t),t))*epsilon^(-eta_order+p+1+zeta_order)-((d*u(x,t)^2-d)*
Gamma(x,t)^2+d-2)*(u(x,t)^2+1)*epsilon^(p+1)*eta*diff(diff(Gamma(x,t),t),x)-((d*u(x,t)^2-d)*Gamma(x,t)^2+d-2)*u(x,t)*epsilon^(p+1)*eta*diff(diff(
Gamma(x,t),t),t)-((d*u(x,t)^2-d)*Gamma(x,t)^2+d-2)*u(x,t)*epsilon^(p+1)*eta*diff(diff(Gamma(x,t),x),x)-u(x,t)*Gamma(x,t)*(d*(u(x,t)^2+3)*Gamma(x,t)^2
+d-2)*epsilon^(p+1)*eta*diff(diff(u(x,t),t),x)-2*Gamma(x,t)*epsilon^(p+1)*eta*(d*u(x,t)^2*Gamma(x,t)^2+d-1)*diff(diff(u(x,t),x),x)-(-1+(u(x,t)^2+1)*
Gamma(x,t)^2)*d*Gamma(x,t)*epsilon^(p+1)*eta*diff(diff(u(x,t),t),t)-2*(((1+(d-1)*u(x,t)^2)*Gamma(x,t)^2+3*d-2)*u(x,t)*Gamma(x,t)^3*diff(u(x,t),x)^2+(
2*((1/2+(d-1/2)*u(x,t)^2)*Gamma(x,t)^2+1/2*d*u(x,t)^2+d-1)*Gamma(x,t)^3*diff(u(x,t),t)+(u(x,t)^2*(u(x,t)-1)*(u(x,t)+1)*(d-1)*Gamma(x,t)^4+((6*d-2)*u(
x,t)^2-d)*Gamma(x,t)^2+2*d-2)*diff(Gamma(x,t),x)+u(x,t)*diff(Gamma(x,t),t)*(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^4+(3/2*d*u(x,t)^2+7/2*d-2)*Gamma(x,t)^2+1
/2*d-1))*diff(u(x,t),x)+d*u(x,t)*Gamma(x,t)^5*diff(u(x,t),t)^2+(u(x,t)*(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^4+(3/2*d*u(x,t)^2+7/2*d-2)*Gamma(x,t)^2+1/2*d
-1)*diff(Gamma(x,t),x)+(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^4+(2*d*u(x,t)^2+3*d-2)*Gamma(x,t)^2-d)*diff(Gamma(x,t),t))*diff(u(x,t),t)+d*Gamma(x,t)*(u(x,t
)-1)*(u(x,t)+1)*(u(x,t)*diff(Gamma(x,t),t)+diff(Gamma(x,t),x))*(u(x,t)*diff(Gamma(x,t),x)+diff(Gamma(x,t),t)))*eta*epsilon^(p+1)+d*(u(x,t)*Gamma(x,t)
^2*(e(x,t)+P(x,t))*diff(u(x,t),x)+Gamma(x,t)^2*(e(x,t)+P(x,t))*diff(u(x,t),t)+u(x,t)*diff(P(x,t),t)+diff(P(x,t),x)))/d = ((1+(u(x,t)^2-1)*Gamma(x,t)^
2)*(u(x,t)*diff(chempot_on_temp(x,t),x)+diff(chempot_on_temp(x,t),t))*temp(x,t)*sigma_Q*epsilon^(q+1/2*m+1/2*abs(m)+1)-Gamma(x,t)*n(x,t)*ee^2*(u(x,t)
-1)*(u(x,t)+1))*diff(phi(x,t),x)/ee, diff(s(x,t),x)*epsilon = (-2*d*sigma_Q*temp(x,t)*Gamma(x,t)*diff(phi(x,t),x)*ee*(u(x,t)*diff(chempot_on_temp(x,t
),t)+diff(chempot_on_temp(x,t),x))*epsilon^(-min(p,q,p+zeta_order-eta_order)+1/2*abs(m)+1/2*m+q)+zeta*d*(u(x,t)*diff(Gamma(x,t),x)+Gamma(x,t)*diff(u(
x,t),x)+diff(Gamma(x,t),t))*((1+(u(x,t)^2-1)*Gamma(x,t)^2)*u(x,t)*diff(Gamma(x,t),x)+(1+(u(x,t)^2-1)*Gamma(x,t)^2)*diff(Gamma(x,t),t)+Gamma(x,t)*((u(
x,t)^2*Gamma(x,t)^2+1)*diff(u(x,t),x)+u(x,t)*Gamma(x,t)^2*diff(u(x,t),t)))*ee^2*epsilon^(-min(p,q,p+zeta_order-eta_order)+p-eta_order+zeta_order)+((
Gamma(x,t)^2-1)*diff(chempot_on_temp(x,t),t)^2+2*u(x,t)*Gamma(x,t)^2*diff(chempot_on_temp(x,t),t)*diff(chempot_on_temp(x,t),x)+diff(chempot_on_temp(x
,t),x)^2*(u(x,t)^2*Gamma(x,t)^2+1))*d*sigma_Q*temp(x,t)^2*epsilon^(-min(p,q,p+zeta_order-eta_order)+abs(m)+m+q)-(-2*eta*((u(x,t)^2*(u(x,t)-1)*(u(x,t)
+1)*(d-1)*Gamma(x,t)^2+(d-1)*u(x,t)^2-1/2*d)*(1+(u(x,t)^2-1)*Gamma(x,t)^2)*diff(Gamma(x,t),x)^2+(2*(1+(u(x,t)^2-1)*Gamma(x,t)^2)*(((d-1)*u(x,t)^2-d+1
)*Gamma(x,t)^2+1/2*d-1)*u(x,t)*diff(Gamma(x,t),t)+2*Gamma(x,t)*(((u(x,t)+1)*(u(x,t)-1)*(1/2+(d-1)*u(x,t)^2)*Gamma(x,t)^4+2*(u(x,t)^2-1/2)*(d-1)*Gamma
(x,t)^2+d-1)*u(x,t)*diff(u(x,t),x)+((u(x,t)+1)*(u(x,t)-1)*(d-1/2)*u(x,t)^2*Gamma(x,t)^4+((3/2*d-1)*u(x,t)^2-1/2*d)*Gamma(x,t)^2+1/2*d)*diff(u(x,t),t)
))*diff(Gamma(x,t),x)+(((d-1)*u(x,t)^2-d+1)*Gamma(x,t)^2+d-1/2*d*u(x,t)^2-1)*(1+(u(x,t)^2-1)*Gamma(x,t)^2)*diff(Gamma(x,t),t)^2+2*Gamma(x,t)*((-1+(u(
x,t)+1)*(u(x,t)-1)*(1/2+(d-1)*u(x,t)^2)*Gamma(x,t)^4+(1+(d-2)*u(x,t)^2)*Gamma(x,t)^2)*diff(u(x,t),x)+((u(x,t)+1)*(u(x,t)-1)*(d-1/2)*Gamma(x,t)^4+(-1/
2*d*u(x,t)^2+3/2*d-1)*Gamma(x,t)^2-1/2*d)*diff(u(x,t),t)*u(x,t))*diff(Gamma(x,t),t)+Gamma(x,t)^2*(((1+(d-1)*u(x,t)^2)*u(x,t)^2*Gamma(x,t)^4+2*(d-1)*u
(x,t)^2*Gamma(x,t)^2+d-1)*diff(u(x,t),x)^2+2*((1/2+(d-1/2)*u(x,t)^2)*Gamma(x,t)^2+d-1)*Gamma(x,t)^2*diff(u(x,t),t)*u(x,t)*diff(u(x,t),x)+d*diff(u(x,t
),t)^2*(-1/2+u(x,t)^2*Gamma(x,t)^4+(-1/2*u(x,t)^2+1/2)*Gamma(x,t)^2)))*epsilon^(-min(p,q,p+zeta_order-eta_order)+p)+d*sigma_Q*Gamma(x,t)^2*diff(phi(x
,t),x)^2*epsilon^(-min(p,q,p+zeta_order-eta_order)+q)*(u(x,t)-1)*(u(x,t)+1))*ee^2)/temp(x,t)/d/ee^2]:

# Test that nondim_orders contains the correct definitions
eval({abs(m),p,q,sigma_Q_order,eta_order,zeta_order}, nondim_orders):
if not type(%,set(numeric)) then
  error("Parameter 'nondim_orders' must be a set/list of the form name=numeric with a name for each of the following: {m,p,q,sigma_Q_order,eta_order,zeta_order}"):
end if:

# Test that abs(m), p, and q are non-negative integers
eval({abs(m),p,q}, nondim_orders):
if not type(%,set(nonnegint)) then
  error("The parameter 'm' in 'nondim_orders' must be a non-zero integer, and the parameters 'p' and 'q' in 'nondim_orders' must by non-negative integers"):
end if:

# Test that m \neq 0
eval(m, nondim_orders):
if % = 0 then
  error("'m' in 'nondim_orders' must be non-zero integer."):
end if:

eval(zeta_order-eta_order, nondim_orders):
if % < 0 then
  error("Equations were derived under the assumption 'zeta_order' \ge 'eta_order'"):
end if:

eval([m, q], nondim_orders):
if %[1] < -1 and %[2] = 0 then
  error("'nondim_orders' with m < -1 and q = 0 is not supported; thermodynamics would require T_1 to depend on n_(1+abs(m)) and P_(1+abs(m)), which breaks the multiple scales approach for abs(m)>1"):
end if:

# Apply nondim_orders to epsilon orders in PDEs
PDEs := simplify(eval(PDEs, nondim_orders)):

cleaner_PDES := [
diff(n(x,t)*Gamma(x,t),t)+diff(n(x,t)*u(x,t)*Gamma(x,t),x) = 0,
diff(Gamma(x,t)^2*(e(x,t)+P(x,t)),t)+diff(Gamma(x,t)^2*u(x,t)*(e(x,t)+P(x,t)),x)-diff(P(x,t),t)=Gamma(x,t)*(-ee*n(x,t))*u(x,t)*E(x,t),
Gamma(x,t)^2*(e(x,t)+P(x,t))*(diff(u(x,t),t)+u(x,t)*diff(u(x,t),x))+(u(x,t)*diff(P(x,t),t)+diff(P(x,t),x))=Gamma(x,t)*(-ee*n(x,t))*(E(x,t)-u(x,t)^2*E(x,t))]:

cleaner_PDES := eval(%,
{E=((x,t)->-diff(phi(x,t),x))}):

# Remove dissipation
eval(PDEs[1..3], {sigma_Q = 0, eta=0, zeta=0}):

expand(% - cleaner_PDES):

CodeTools:-Test(%, [0,0,0], label="Verify form of PDEs", quiet=true):

# Insert electric potential
phi = unapply(AA/(-ee)*(Gamma(x,t)*n(x,t) + 1/3*epsilon*d1*d2*diff(Gamma(x,t)*n(x,t),x,x)),x,t):
PDEs := eval(PDEs,%):

# Insert definition of Gamma
Gamma = unapply(1/sqrt(1-u(x,t)^2),x,t):
PDEs := eval(PDEs,%):

# Insert definition of chempot_on_temp
chempot_on_temp = unapply(chempot(x,t)/temp(x,t),x,t):
PDEs := eval(PDEs, %):

# AA:=4*Pi*ee^2*d1*d2/kappa/(d1+d2):

# Choose units where ee=1 (we already set hbar=k__B=v__F=1); don't do this
# earlier since it's useful (in top "calculate equations") section to see where
# ee is, especially since it can sometime be unclear if ee is postive or
# negative.
PDEs := eval(PDEs,ee=1):

ind_vars := {x}:

pdes, IndDefs, DepDefs, new_ind_vars := MultipleScales:-GenEqns(PDEs,
  e_order, ':-scales'=t, ':-scales_order'=t_num,
  ':-leading_order_eqs'=leading_order_eqs, ':-ind_vars'=ind_vars)[1..4]:

return pdes:

end proc:

end module:

local Latex := proc(eq::equation)
  convert(eq,diff):
  latex(%):
  printf(`\\\\\n`):
end proc:

#Graphene:-UnitTests(quiet_tests=true):