diff --git a/README.md b/README.md
index 536fe98..addfaa3 100644
--- a/README.md
+++ b/README.md
@@ -1,9 +1,9 @@
-:green_book: The [manual](https://github.com/Leedehai/typst-physics/blob/master/physica-manual.pdf).
+:green_book: The [manual](https://github.com/Leedehai/typst-physics/blob/v0.9.1/physica-manual.pdf).
 <p align="center">
 <img width="545" alt="logo" src="https://github.com/Leedehai/typst-physics/assets/18319900/ed86198a-8ddb-4473-aed3-8111d5ecde60">
 </p>
 
-# The physica package for Typst
+# The physica package for Typst (v0.9.1)
 
 [![CI](https://github.com/Leedehai/typst-physics/actions/workflows/ci.yml/badge.svg)](https://github.com/Leedehai/typst-physics/actions/workflows/ci.yml)
 [![Latest release](https://img.shields.io/github/v/release/Leedehai/typst-physics.svg?color=gold)][latest-release]
@@ -27,7 +27,7 @@ natural sciences, including:
 
 ## A quick look
 
-See the [manual](https://github.com/Leedehai/typst-physics/blob/master/physica-manual.pdf) for more details and examples.
+See the [manual](https://github.com/Leedehai/typst-physics/blob/v0.9.1/physica-manual.pdf)for more details and examples.
 
 ![demo-quick](https://github.com/Leedehai/typst-physics/assets/18319900/4a9f40df-f753-4324-8114-c682d270e9c7)
 
@@ -89,7 +89,7 @@ typst 0.10.0 (70ca0d25)
 
 ## Manual
 
-See the manual [physica-manual.pdf](https://github.com/Leedehai/typst-physics/blob/master/physica-manual.pdf) for a more comprehensive coverage, a PDF file
+See the [manual](https://github.com/Leedehai/typst-physics/blob/v0.9.1/physica-manual.pdf) for a more comprehensive coverage, a PDF file
 generated directly with the [Typst](https://typst.app) binary.
 
 To regenerate the manual, use command
@@ -107,8 +107,9 @@ request. If it is large, the best first step is creating an issue and let us
 explore the design together. Some features might warrant a package on its own.
 
 * Testing: currently testing is done by closely inspecting the generated
-[physica-manual.pdf](https://github.com/Leedehai/typst-physics/blob/master/physica-manual.pdf). This does not scale well. I plan to add programmatic
-testing by comparing rendered pictures with golden images.
+[manual](https://github.com/Leedehai/typst-physics/blob/v0.9.1/physica-manual.pdf).
+This does not scale well. I plan to add programmatic testing by comparing
+rendered pictures with golden images.
 
 ## License
 
diff --git a/changelog.md b/changelog.md
new file mode 100644
index 0000000..8178e9e
--- /dev/null
+++ b/changelog.md
@@ -0,0 +1,66 @@
+# Changelog
+
+## 0.9.2 (pending)
+
+* Let `braket()` take an optional third argument to render `braket(u,A,v)` as
+`⟨u|A|v⟩`. The one argument case `braket(a)` is still rendered as `⟨a|a⟩` and
+the two argument case `braket(u,v)` is still rendered as `⟨u|v⟩`.
+* Let `expval()` take an optional second argument to render `expval(A,a)` as
+`⟨a|A|a⟩`. The one argument case `expval(A)` is still rendered as `⟨A⟩`.
+* **(breaking)** Differentiated big-O `Order(...)` and little-o `order(...)`.
+* Let Taylor series term `taylorterm(...)` automatically add parenthesis, so
+that e.g. `tarlorterm(f,x,1+a,n-1)`'s `1+a` and `n-1` can be put inside
+parenthesis when needed.
+* Fixed a rendering issue of `Set(...)` when it contains tall contents.
+* Removed the `box(..)` layer from `tensor(...)`'s phantom index.
+* Added 2D and 3D rotation matrix `rot2mat(...)` and `rot3mat(...)`.
+* Added Gram matrix `grammat(...)`.
+* Removed `gradient` and `divergence`, since most users will use the abbreviated
+names `grad` and `div`. This prevents name collisions when users do wildcard
+importing, especially with Typst's built-in `gradient` that shows color
+gradients on texts.
+
+## 0.9.1
+
+* Added show rules `super-T-as-transpose` and `super-plus-as-dagger`, so that
+`..^T` and `..^+` can be rendered as a transpose operator and a dagger
+(i.e. conjugate transpose) operator, respectively, just like handwriting.
+* Adjusted spacing for `dd()`, following advice in TeXBook Chapter 18. The
+intra spacing can be disabled by a new optional argument `compact:#true`.
+* Fixed a spacing issue of `vecrow(..)`.
+* **(breaking)** Let `iprod(..)` be rendered in the more familiar inner
+product `⟨a,b⟩`, instead of `⟨a|b⟩`--that can be done with `braket(a,b)`.
+* Added optional argument `big:#true` to Jacobian matrix and Hessian matrix, so
+that the fractions can be rendered in an "uncramped" form. The default is the
+"cramped" form just like LaTeX.
+* Added `taylorterm(...)` to display terms in the Taylor series.
+
+## 0.9.0
+
+Changed the minimum Typst version to 0.10.0.
+
+* Fixed the appearance of `braket`, `ketra`, and `mel`.
+* Fixed Hessian matrix `hmax(...)`.
+
+## 0.8.1
+
+## 0.8.0
+
+* Added CI.
+* **(breaking)** Let `va(...)` (vector arrow) be not bold, according to the ISO
+standard.
+
+## 0.7.5
+
+* First appeared in the official package collection
+[typst/packages](https://github.com/typst/packages).
+
+## 0.7.4
+
+* Renamed from `physics` to `physica`, meaning _natural sciences_.
+
+## 0.7.1
+
+## 0.6
+
+## 0.5
diff --git a/physica-manual.pdf b/physica-manual.pdf
index f062981..dd13c91 100644
Binary files a/physica-manual.pdf and b/physica-manual.pdf differ
diff --git a/physica-manual.typ b/physica-manual.typ
index 31a584d..0cb3e61 100644
--- a/physica-manual.typ
+++ b/physica-manual.typ
@@ -109,10 +109,15 @@ All symbols need to be used in *math mode* `$...$`.
   [`norm(phi(x))`   #sym.arrow $norm(phi(x))$],
   [norm],
 
+  [`Order(`_content_`)`],
+  [],
+  [`Order(x^2)`   #sym.arrow $Order(x^2)$],
+  [big O],
+
   [`order(`_content_`)`],
   [],
-  [`order(x^2)`   #sym.arrow $order(x^2)$],
-  [order of magnitude],
+  [`order(1)`   #sym.arrow $order(1)$],
+  [small O],
 
   [`Set(`_content_`)`],
   [],
@@ -134,9 +139,9 @@ All symbols need to be used in *math mode* `$...$`.
   [`expval`],
   [
     `expval(u)`  #sym.arrow $expval(u)$ \
-    `expval(f/N)`  #sym.arrow $expval(f/N)$ \
+    `expval(p,psi)`  #sym.arrow $expval(p,psi)$ \
   ],
-  [expectation value],
+  [expectation value, also see bra-ket @dirac-braket below],
 )
 
 == Vector notations
@@ -182,13 +187,13 @@ All symbols need to be used in *math mode* `$...$`.
   [`va(a),va(mu_1)` #sym.arrow $va(a),va(mu_1)$],
   [vector, arrow \ #sub[(not bold: see ISO 80000-2:2019)]],
 
-  [`gradient`],
   [`grad`],
+  [],
   [`grad f` #sym.arrow $grad f$],
   [gradient],
 
-  [`divergence`],
   [`div`],
+  [],
   [`div vb(E)` #sym.arrow $div vb(E)$],
   [divergence],
 
@@ -225,6 +230,8 @@ All symbols need to be used in *math mode* `$...$`.
 
 #v(1em)
 
+=== Determinant, (anti-)diagonal, identity, zero matrix
+
 #table(
   columns: (auto, auto, auto, auto),
   align: left,
@@ -279,24 +286,11 @@ All symbols need to be used in *math mode* `$...$`.
     #hl(`zmat(3,delim:"[")`) #sym.arrow \ $zmat(3,delim:"[")$
   ],
   [zero matrix],
-
-  [`jacobianmatrix(`...`)`],
-  [`jmat`],
-  [See below],
-  [Jacobian matrix],
-
-  [`hessianmatrix(`...`)`],
-  [`hmat`],
-  [See below],
-  [Hessian matrix],
-
-  [`xmatrix(`...`)`],
-  [`xmat`],
-  [See below],
-  [Matrix built with an element building function],
 )
 
-Jacobian matrix: `jacobianmatrix(`...`)`, i.e. `jmat(`...`)`.
+=== Jacobian matrix
+
+`jacobianmatrix(`...`)`, i.e. `jmat(`...`)`.
 
 #table(
   columns: (25%, auto, auto),
@@ -328,7 +322,9 @@ Jacobian matrix: `jacobianmatrix(`...`)`, i.e. `jmat(`...`)`.
   ],
 )
 
-Hessian matrix: `hessianmatrix(`...`)`, i.e. `hmat(`...`)`.
+=== Hessian matrix
+
+`hessianmatrix(`...`)`, i.e. `hmat(`...`)`.
 
 #table(
   columns: (25%, auto, auto),
@@ -360,7 +356,10 @@ Hessian matrix: `hessianmatrix(`...`)`, i.e. `hmat(`...`)`.
   ],
 )
 
-Matrix built with an element building function: `xmatrix(`_m, n, func_`)`, i.e. `xmat(`...`)`. The element building function _func_ takes two integers which are the row and column numbers starting from 1.
+=== Matrix with an element builder
+
+`xmatrix(`_m, n, func_`)`, i.e. `xmat(`...`)`. The element building function
+_func_ takes two integers which are the row and column numbers starting from 1.
 
 #table(
   columns: (auto, auto),
@@ -385,7 +384,64 @@ xmat(2, 2, #g)`)
   ],
 )
 
-== Dirac braket notations
+=== Rotation matrices, 2D and 3D
+
+#table(
+  columns: (auto, auto, auto),
+  align: center,
+  stroke: none,
+  column-gutter: 1em,
+
+  [
+    #hl(`rot2mat(theta)`)
+    $ rot2mat(theta) $
+  ],
+  [
+    #hl(`rot2mat(-a/2,delim:"[")`)
+    $ rot2mat(-a/2, delim:"[") $
+  ],
+  [
+    #hl(`rot2mat(display(a/2),delim:"[")`)
+    $ rot2mat(display(a/2),delim:"[") $
+  ],
+
+  [
+    #hl(`rot3xmat(theta)`)
+    $ rot3xmat(theta) $
+  ],
+  [
+    #hl(`rot3ymat(45^degree)`)
+    $ rot3ymat(45^degree) $
+  ],
+  [
+    #hl(`rot3zmat(theta,delim:"[")`)
+    $ rot3zmat(theta,delim:"[") $
+  ],
+)
+
+=== Gram matrix
+
+#table(
+  columns: (auto, auto, auto),
+  align: center,
+  stroke: none,
+  column-gutter: 1em,
+
+  [
+    #hl(`grammat(alpha,beta)`)
+    $ grammat(alpha, beta) $
+  ],
+  [
+    #hl(`grammat(v_1,v_2,v_3, delim:"[")`)
+    $ grammat(v_1,v_2,v_3, delim:"[") $
+  ],
+  [
+    #hl(`grammat(v_1,v_2, norm:#true)`)
+    $ grammat(v_1,v_2, norm:#true) $
+  ],
+)
+
+== Dirac braket notations <dirac-braket>
 
 #v(1em)
 
@@ -411,37 +467,36 @@ xmat(2, 2, #g)`)
   ],
   [ket],
 
-  [`expval(`_content_`)`],
+  [`braket(`..`)`],
   [],
   [
-    `expval(u)`  #sym.arrow $expval(u)$ \
-    `expval(vec(1,2))` #sym.arrow $expval(vec(1,2))$
+    `braket(a), braket(u, v)` \ #sym.arrow $braket(a), braket(u, v)$ \
+    `braket(psi,A/N,phi)` #sym.arrow $braket(psi,A/N,phi)$
   ],
-  [expectation],
+  [braket, with 1, 2, or 3 arguments],
 
-  [`braket(`_a_, _b_`)`],
+  [`ketbra(`..`)`],
   [],
   [
-    `braket(a), braket(u, v)` \ #sym.arrow $braket(a), braket(u, v)$ \
-    `braket(vec(1,2), b)` #sym.arrow $braket(vec(1,2), b)$
+    `ketbra(a), ketbra(u, v)` \ #sym.arrow $ketbra(a), ketbra(u, v)$ \
+    `ketbra(a/N, b)` #sym.arrow $ketbra(a/N, b)$
   ],
-  [braket],
+  [ketbra, with 1 or 2 arguments],
 
-  [`ketbra(`_a_, _b_`)`],
+  [`expval(`_content_`)`],
   [],
   [
-    `ketbra(a), ketbra(u, v)` \ #sym.arrow $ketbra(a), ketbra(u, v)$ \
-    `ketbra(vec(1,2), b)` #sym.arrow $ketbra(vec(1,2), b)$
+    `expval(u)`  #sym.arrow $expval(u)$ \
+    `expval(A,psi)` #sym.arrow $expval(A,psi)$
   ],
-  [ketbra],
+  [expectation],
 
-  [`matrixelement(`_n_, _M_, _m_`)`],
+  [`matrixelement(`..`)`],
   [`mel`],
   [
-    `mel(n, diff_nu H, m)` \ #sym.arrow $mel(n, diff_nu H, m)$ \
-    `mel(vec(U,V),A,m)` #sym.arrow $mel(vec(U,V),A,m)$
+    `mel(n, diff_nu H, m)` \ #sym.arrow $mel(n, diff_nu H, m)$
   ],
-  [matrix element],
+  [matrix element, same as `braket(n,M,n)`],
 )
 
 == Math functions
@@ -962,8 +1017,8 @@ Function: `tensor(`_symbol_, \*_args_`)`.
     $ tensor(AA,+a,+b,-c,-d,+e,-f,+g,-h) $
   ],
   [
-    *(7)* #hl(`tensor(R, -a, -b, -c, +d)`) \
-    $ tensor(R, -a, -b, -c, +d) $
+    *(7)* #hl(`tensor(R, -a, -b, +d)`) \
+    $ tensor(R, -a, -b, +d) $
   ],
   [
     *(8)* #hl(`tensor(T,+1,-I(1,-1),+a_bot,-+,+-)`) \
@@ -1001,11 +1056,11 @@ Function: `isotope(`_element_, _a_: ..., _z_: ...`)`.
   ],
 )
 
-*(3)* #hl(`isotope("Bi",a:211,z:83) --> isotope("Tl",a:207,z:81) + isotope("He",a:4,z:2)`)
-$ isotope("Bi",a:211,z:83) --> isotope("Tl",a:207,z:81) + isotope("He",a:4,z:2) $
+*(3)* #hl(`isotope("Bi",a:211,z:83) -> isotope("Tl",a:207,z:81) + isotope("He",a:4,z:2)`)
+$ isotope("Bi",a:211,z:83) -> isotope("Tl",a:207,z:81) + isotope("He",a:4,z:2) $
 
-*(4)* #hl(`isotope("Tl",a:207,z:81) --> isotope("Pb",a:207,z:82) + isotope(e,a:0,z:-1)`)
-$ isotope("Tl",a:207,z:81) --> isotope("Pb",a:207,z:82) + isotope(e,a:0,z:-1) $
+*(4)* #hl(`isotope("Tl",a:207,z:81) -> isotope("Pb",a:207,z:82) + isotope(e,a:0,z:-1)`)
+$ isotope("Tl",a:207,z:81) -> isotope("Pb",a:207,z:82) + isotope(e,a:0,z:-1) $
 
 === The n-th term in Taylor series
 
@@ -1016,8 +1071,8 @@ Function: `taylorterm(`_func_, _x_, _x0_, _idx_`)`.
 - _x_: the variable name e.g. `x`,
 - _x0_: the variable value at the expansion point e.g. `x_0`, `(1+a)`,
 - _idx_: the index of the term, e.g. `0`, `1`, `2`, `n`, `(n+1)`.
-If _x0_ or _idx_ is in a parenthesis e.g. `(1+k)`, then the function
-automatically removes the outer parenthesis where appropriate.
+If _x0_ or _idx_ is an add/sub sequence e.g. `-a`, `a+b`, then the function
+automatically adds a parenthesis where appropriate.
 
 #align(center, [*Examples*])
 
@@ -1034,20 +1089,20 @@ automatically removes the outer parenthesis where appropriate.
     $ taylorterm(f,x,x_0,1) $
   ],
   [
-    *(3)* #hl(`taylorterm(f,x,(1+a),2)`) \
-    $ taylorterm(f,x,(1+a),2) $
+    *(3)* #hl(`taylorterm(F,x^nu,x^nu_0,n)`) \
+    $ taylorterm(F,x^nu,x^nu_0,n) $
   ],
   [
     *(4)* #hl(`taylorterm(f,x,x_0,n)`) \
     $ taylorterm(f,x,x_0,n) $
   ],
   [
-    *(5)* #hl(`taylorterm(F,x^nu,x^nu_0,n)`) \
-    $ taylorterm(F,x^nu,x^nu_0,n) $
+    *(5)* #hl(`taylorterm(f,x,1+a,2)`) \
+    $ taylorterm(f,x,1+a,2) $
   ],
   [
-    *(6)* #hl(`taylorterm(f_p,x,x_0,(n+1))`) \
-    $ taylorterm(f_p,x,x_0,(n+1)) $
+    *(6)* #hl(`taylorterm(f_p,x,x_0,n-1)`) \
+    $ taylorterm(f_p,x,x_0,n-1) $
   ],
 )
 
diff --git a/physica.typ b/physica.typ
index 4b78249..7641823 100644
--- a/physica.typ
+++ b/physica.typ
@@ -3,8 +3,32 @@
 // Repository: https://github.com/Leedehai/typst-physics
 // Please see physica-manual.pdf for user docs.
 
+// Returns whether a Content object is an add/sub sequence, e.g. -a, a+b, a-b.
+// The caller is responsible for ensuring the input is a Content object.
+#let __is_add_sub_sequence(content) = {
+  if not content.has("children") { return false }
+
+  let impl(seq) = {
+    // Only check the top level, don't descend into the child, since we don't
+    // care if the child is a parenthesis group that contains +/-.
+    for child in seq.at("children") {
+      if child == [+] or child == [#sym.minus] { return true }
+    }
+    return false
+  }
+
+  // Before my https://github.com/typst/typst/pull/3063:
+  return impl(content)
+  // After my https://github.com/typst/typst/pull/3063:
+  // if content.func() == math.math-style {
+  //   return impl(content.at("body"))
+  // } else {
+  //   return impl(content)
+  // }
+}
+
 // Returns whether a Content object holds an integer. The caller is responsible
-// for ensuring the input argument is a Content object.
+// for ensuring the input is a Content object.
 #let __content_holds_number(content) = {
   return content.func() == text and regex("^\d+$") in content.text
 }
@@ -123,28 +147,37 @@
 
 // == Braces
 
-#let Set(..sink) = style(styles => {
+#let Set(..sink) = {
   let args = sink.pos()  // array
   let expr = args.at(0, default: none)
   let cond = args.at(1, default: none)
-  let height = measure($ expr cond $, styles).height;
-  let phantom = box(height: height, width: 0pt, inset: 0pt, stroke: none);
 
   if expr == none {
-    if cond == none { ${}$ } else { ${lr(|phantom#h(0pt))#cond}$ }
+    if cond == none { ${}$ } else { ${mid(|) #cond}$ }
   } else {
-    if cond == none { ${#expr}$ } else { ${#expr lr(|phantom#h(0pt))#cond}$ }
+    if cond == none { ${#expr}$ } else { ${#expr mid(|) #cond}$ }
   }
-})
+}
 
-#let order(content) = $cal(O)(content)$
+#let Order(content) = $cal(O)(content)$
+#let order(content) = $cal(o)(content)$
 
 #let evaluated(content) = {
   $lr(zwj#content|)$
 }
 #let eval = evaluated
 
-#let expectationvalue(f) = $lr(angle.l #f angle.r)$
+#let expectationvalue(..sink) = {
+  let args = sink.pos()  // array
+  let expr = args.at(0, default: none)
+  let func = args.at(1, default: none)
+
+  if func == none {
+    $lr(angle.l expr angle.r)$
+  } else {
+    $lr(angle.l func#h(0pt)mid(|)#h(0pt)expr#h(0pt)mid(|)#h(0pt)func angle.r)$
+  }
+}
 #let expval = expectationvalue
 
 // == Vector notations
@@ -196,14 +229,9 @@
 #let vectorarrow(a) = __vector(a, math.arrow, false)
 #let va = vectorarrow
 
-#let gradient = $bold(nabla)$
-#let grad = gradient
-
-#let divergence = $bold(nabla)dot.c$
-#let div = divergence
-
+#let grad = $bold(nabla)$
+#let div = $bold(nabla)dot.c$
 #let curl = $bold(nabla)times$
-
 #let laplacian = $nabla^2$
 
 #let dotproduct = $dot$
@@ -273,11 +301,12 @@
 #let admat = antidiagonalmatrix
 
 #let identitymatrix(order, delim:"(", fill:none) = {
-  let order_num = 1
-  if type(order) == content and __content_holds_number(order) {
-    order_num = int(order.text)
+  let order_num = if type(order) == content and __content_holds_number(order) {
+    int(order.text)
+  } else if type(order) == "integer" {
+    order
   } else {
-    panic("the order shall be an integer, e.g. 2")
+    panic("imat/identitymatrix: the order shall be an integer, e.g. 2")
   }
 
   let ones = range(order_num).map((i) => 1)
@@ -286,11 +315,12 @@
 #let imat = identitymatrix
 
 #let zeromatrix(order, delim:"(") = {
-  let order_num = 1
-  if type(order) == content and __content_holds_number(order) {
-    order_num = int(order.text)
+  let order_num = if type(order) == content and __content_holds_number(order) {
+    int(order.text)
+  } else if type(order) == "integer" {
+    order
   } else {
-    panic("the order shall be an integer, e.g. 2")
+    panic("zmat/zeromatrix: the order shall be an integer, e.g. 2")
   }
 
   let ones = range(order_num).map((i) => 0)
@@ -333,18 +363,22 @@
 #let hmat = hessianmatrix
 
 #let xmatrix(m, n, func, delim:"(") = {
-  let rows = none
-  if type(m) == content and __content_holds_number(m) {
-    rows = int(m.text)
+  let rows = if type(m) == content and __content_holds_number(m) {
+    int(m.text)
+  } else if type(m) == "integer" {
+    m
   } else {
-    panic("the first argument shall be an integer, e.g. 2")
+    panic("xmat/xmatrix: the first argument shall be an integer, e.g. 2")
   }
-  let cols = none
-  if type(n) == content and __content_holds_number(m) {
-    cols = int(n.text)
+
+  let cols = if type(n) == content and __content_holds_number(m) {
+    int(n.text)
+  } else if type(n) == "integer" {
+    n
   } else {
-    panic("the second argument shall be an integer, e.g. 2")
+    panic("xmat/xmatrix: the second argument shall be an integer, e.g. 2")
   }
+
   assert(
     type(func) == function,
     message: "func shall be a function (did you forget to add a preceding '#' before the function name)?"
@@ -361,6 +395,56 @@
 }
 #let xmat = xmatrix
 
+#let rot2mat(theta, delim:"(") = {
+  let operand = if type(theta) == "content" and __is_add_sub_sequence(theta) {
+    $(theta)$
+  } else { theta }
+  $mat(cos operand, -sin operand;
+       sin operand, cos operand; delim: delim)$
+}
+
+#let rot3xmat(theta, delim:"(") = {
+  let operand = if type(theta) == "content" and __is_add_sub_sequence(theta) {
+    $(theta)$
+  } else { theta }
+  $mat(1, 0,           0;
+       0, cos operand, -sin operand;
+       0, sin operand, cos operand; delim: delim)$
+}
+
+#let rot3ymat(theta, delim:"(") = {
+  let operand = if type(theta) == "content" and __is_add_sub_sequence(theta) {
+    $(theta)$
+  } else { theta }
+  $mat(cos operand,  0, sin operand;
+       0,            1, 0;
+       -sin operand, 0, cos operand; delim: delim)$
+}
+
+#let rot3zmat(theta, delim:"(") = {
+  let operand = if type(theta) == "content" and __is_add_sub_sequence(theta) {
+    $(theta)$
+  } else { theta }
+  $mat(cos operand, -sin operand, 0;
+       sin operand, cos operand,  0;
+       0,           0,            1; delim: delim)$
+}
+
+#let grammat(..sink) = {
+  let vs = sink.pos()  // array
+  let delim = sink.named().at("delim", default: "(")
+  let asnorm = sink.named().at("norm", default: false)
+
+  xmat(vs.len(), vs.len(), (i,j) => {
+    if (i == j and (not asnorm)) or i != j {
+      iprod(vs.at(i - 1), vs.at(j - 1))
+    } else {
+      let v = vs.at(i - 1)
+      $norm(#v)^2$
+    }
+  }, delim: delim)
+}
+
 // == Dirac braket notations
 
 #let bra(f) = $lr(angle.l #f|)$
@@ -368,12 +452,16 @@
 
 #let braket(..sink) = style(styles => {
   let args = sink.pos()  // array
-  assert(args.len() == 1 or args.len() == 2, message: "expecting 1 or 2 args")
 
-  let bra = args.at(0)
-  let ket = args.at(1, default: bra)
+  let bra = args.at(0, default: none)
+  let ket = args.at(-1, default: bra)
 
-  $ lr(angle.l bra#h(0pt)mid(bar.v)#h(0pt)ket angle.r) $
+  if args.len() <= 2 {
+    $ lr(angle.l bra#h(0pt)mid(|)#h(0pt)ket angle.r) $
+  } else {
+    let middle = args.at(1)
+    $ lr(angle.l bra#h(0pt)mid(|)#h(0pt)middle#h(0pt)mid(|)#h(0pt)ket angle.r) $
+  }
 })
 
 #let ketbra(..sink) = style(styles => {
@@ -383,11 +471,11 @@
   let ket = args.at(0)
   let bra = args.at(1, default: ket)
 
-  $ lr(bar.v ket#h(0pt)mid(angle.r#h(0pt)angle.l)#h(0pt)bra bar.v) $
+  $ lr(|ket#h(0pt)mid(angle.r#h(0pt)angle.l)#h(0pt)bra|) $
 })
 
 #let matrixelement(n, M, m) = style(styles => {
-  $ lr(angle.l #n#h(0pt)mid(bar.v)#h(0pt)#M#h(0pt)mid(bar.v)#h(0pt)#m angle.r) $
+  $ lr(angle.l #n#h(0pt)mid(|)#h(0pt)#M#h(0pt)mid(|)#h(0pt)#m angle.r) $
 })
 
 #let mel = matrixelement
@@ -713,7 +801,7 @@
   let args = sink.pos()
 
   let (uppers, lowers) = ((), ())  // array, array
-  let hphantom(s) = { hide(box(height: 0em, s)) }  // Like Latex's \hphantom
+  let hphantom(s) = { hide($#s$) }  // Like Latex's \hphantom
 
   for i in range(args.len()) {
     let arg = args.at(i)
@@ -751,21 +839,17 @@
 }
 
 #let taylorterm(fn, xv, x0, idx) = {
-  let noparen(expr) = {
-    if type(expr) == content and expr.func() == math.lr {
-      let children = expr.at("body").at("children")
-      children.slice(1, children.len() - 1).join()
-    } else {
-      expr
-    }
+  let maybeparen(expr) = {
+    if __is_add_sub_sequence(expr) { $(expr)$ }
+    else { expr }
   }
 
   if idx == [0] or idx == 0 {
-    $fn (noparen(x0))$
+    $fn (x0)$
   } else if idx == [1] or idx == 1 {
-    $fn^((1)) (noparen(x0))(xv - x0)$
+    $fn^((1)) (x0)(xv - maybeparen(x0))$
   } else {
-    $frac(fn^((noparen(idx))) (noparen(x0)), idx !)(xv - x0)^noparen(idx)$
+    $frac(fn^((idx)) (x0), maybeparen(idx) !)(xv - maybeparen(x0))^idx$
   }
 }