diff --git a/del-geo-core/src/quaternion.rs b/del-geo-core/src/quaternion.rs
index ef0b6e7..8f9260c 100644
--- a/del-geo-core/src/quaternion.rs
+++ b/del-geo-core/src/quaternion.rs
@@ -107,11 +107,19 @@ where
     ]
 }
 
-pub fn from_axisangle(a: &[f32; 3]) -> [f32; 4] {
-    let half = 0.5;
+pub fn from_axisangle<Real>(a: &[Real; 3]) -> [Real; 4]
+where
+    Real: num_traits::Float,
+{
+    let half = Real::from(0.5).unwrap();
     let sqlen = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
-    if sqlen < 1.0e-10 {
-        return [half * a[0], half * a[1], half * a[2], 1. - 0.125 * sqlen];
+    if sqlen < Real::epsilon() {
+        return [
+            half * a[0],
+            half * a[1],
+            half * a[2],
+            -sqlen * Real::from(0.125).unwrap() + Real::one(),
+        ];
     }
     let lena = sqlen.sqrt();
     [
@@ -127,9 +135,12 @@ pub fn identity() -> [f64; 4] {
 }
 
 /// return rotation around axis with radian
-pub fn around_axis(a: &[f64; 3], rad: f64) -> [f64; 4] {
+pub fn around_axis<Real>(a: &[Real; 3], rad: Real) -> [Real; 4]
+where
+    Real: num_traits::Float,
+{
     let v = crate::vec3::normalized(a);
-    let half = rad * 0.5f64;
+    let half = rad / Real::from(2).unwrap();
     let sin = half.sin();
     [v[0] * sin, v[1] * sin, v[2] * sin, half.cos()]
 }
diff --git a/del-geo-core/src/spherical_harmonics.rs b/del-geo-core/src/spherical_harmonics.rs
index 754d8b1..3f748b8 100644
--- a/del-geo-core/src/spherical_harmonics.rs
+++ b/del-geo-core/src/spherical_harmonics.rs
@@ -264,17 +264,11 @@ pub fn sph_coeff_buffer(n: i8, x: f64, y: f64, z: f64) -> [f64; 100] {
 /// Calculate the coefficients of Legendre Polynomials in all orders <= n.
 /// Try to access the coefficient of x^m in P_l(x) by result[l][m].
 fn legendre_coeff_vec(n: u64) -> Vec<Vec<f64>> {
-    let mut res = Vec::new();
+    let mut res = vec![vec![1.0], vec![0.0, 1.0]];
 
-    let p0 = vec![1.0];
-    res.push(p0);
-
-    let p1 = vec![0.0, 1.0];
-    res.push(p1);
-
-    for i in 2..n + 1 {
-        let mut p: Vec<f64> = Vec::new();
-        for j in 0..i + 1 {
+    for i in 2..=n {
+        let mut p = Vec::new();
+        for j in 0..=i {
             if (i - j) % 2 != 0 {
                 p.push(0.0);
             } else if j == 0 {
@@ -301,15 +295,14 @@ fn legendre_coeff_vec(n: u64) -> Vec<Vec<f64>> {
 
 /// Calculate the factorial of a number.
 /// Here needs a lot of optimization to avoid overflow. Now the maximum parameter is 11, which is enough for the current use.
-fn factorial(n: u128) -> u128 {
+const fn factorial(n: u128) -> u128 {
+    assert!(n <= 34, "The factorial is too large to calculate.");
+    // use recursion because const_for is unstable https://github.com/rust-lang/rust/issues/87575
     if n == 0 {
-        return 1;
-    }
-    let mut res = 1;
-    for i in 1..n + 1 {
-        res *= i;
+        1
+    } else {
+        n * factorial(n - 1)
     }
-    res
 }
 
 /// Calculate the coefficient of the term in the Legendre Polynomial.
@@ -320,10 +313,10 @@ pub fn get_legendre_poly_term_coeff(func_order: u32, term_order: u32) -> f64 {
     } else {
         let k = (func_order - term_order) / 2;
         let mol = if k % 2 != 0 { -1.0 } else { 1.0 }
-            * factorial((func_order + term_order) as u128) as f64;
+            * factorial((func_order + term_order).into()) as f64;
         let den = (2_u32.pow(func_order) as u128
-            * factorial(k as u128)
-            * factorial(term_order as u128)
+            * factorial(k.into())
+            * factorial(term_order.into())
             * factorial((k + term_order).into())) as f64;
         mol / den
     }
@@ -337,24 +330,25 @@ fn calculate_assoc_legendre_poly(l: u64, m: i64, x: f64) -> f64 {
     let sign = m.signum();
     let legendre = legendre_coeff_vec(l);
     let legendre_coeff_vec = &legendre[l as usize];
-    let mut sum_coeff: Vec<f64> = Vec::new();
-    for i in m_abs as usize..l as usize + 1 {
+    let mut sum_coeff = Vec::new();
+    for i in m_abs as usize..=l as usize {
         let mut tmp = legendre_coeff_vec[i];
         tmp *= factorial(i as u128) as f64 / factorial((i as u32 - m_abs as u32).into()) as f64;
         sum_coeff.push(tmp);
     }
-    let mut sum = 0.0;
-    for i in 0..sum_coeff.len() {
-        sum += sum_coeff[i] * f64::powi(x, i as i32);
-    }
+    let sum: f64 = sum_coeff
+        .iter()
+        .enumerate()
+        .map(|(i, &coeff)| coeff * x.powi(i as i32))
+        .sum();
     let mut res =
         if m_abs % 2 == 0 { 1.0 } else { -1.0 } * (1.0 - x * x).powf(m_abs as f64 / 2.0) * sum;
     res *= if sign == 1 {
         1.0
     } else {
         (if m.abs() % 2 == 0 { 1.0 } else { -1.0 })
-            * factorial((l - m.unsigned_abs()) as u128) as f64
-            / factorial((l + m.unsigned_abs()) as u128) as f64
+            * factorial((l - m.unsigned_abs()).into()) as f64
+            / factorial((l + m.unsigned_abs()).into()) as f64
     } as f64;
     res
 }
diff --git a/del-geo-core/src/view_rotation.rs b/del-geo-core/src/view_rotation.rs
index dd38caf..a12feb0 100644
--- a/del-geo-core/src/view_rotation.rs
+++ b/del-geo-core/src/view_rotation.rs
@@ -16,7 +16,7 @@ impl Trackball {
     pub fn camera_rotation(&mut self, cursor_dx: f64, cursor_dy: f64) {
         let dx = cursor_dx as f32;
         let dy = cursor_dy as f32;
-        let a: f32 = (dx * dx + dy * dy).sqrt();
+        let a = (dx * dx + dy * dy).sqrt();
         if a == 0.0 {
             return;
         }