From 167d41be0a646b438aed51ab90ae76efe33607cf Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 22 Apr 2024 10:58:11 +0200
Subject: [PATCH 01/10] Add elevation rasters to scene attributes #5

---
 src/index.ts | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/index.ts b/src/index.ts
index 6df7991..5118fdf 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -20,6 +20,7 @@ import { rayTracingWebGL } from './rayTracingWebGL.js';
 export default class Scene {
   simulationGeometries: Array<BufferGeometry>;
   shadingGeometries: Array<BufferGeometry>;
+  elevationRasters: Array<Array<number[]>>;
   latitude: number;
   longitude: number;
 
@@ -31,6 +32,7 @@ export default class Scene {
   constructor(latitude: number, longitude: number) {
     this.simulationGeometries = [];
     this.shadingGeometries = [];
+    this.elevationRasters = [];
     this.latitude = latitude;
     this.longitude = longitude;
   }
@@ -57,6 +59,10 @@ export default class Scene {
     this.shadingGeometries.push(geometry);
   }
 
+  addElevationRaster(raster: number[][]) {
+    this.elevationRasters.push(raster);
+  }
+
   /** @ignore */
   refineMesh(mesh: BufferGeometry, maxLength: number): BufferGeometry {
     const positions = mesh.attributes.position.array.slice();

From 9d8d690a32df7098e6c3c34779506db0353835b7 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 22 Apr 2024 13:14:31 +0200
Subject: [PATCH 02/10] Create function to calculate elevationAngles #5

Based on a 2D grid which will be given as GeoTIFFs on a running website.
---
 src/elevation.ts        | 53 +++++++++++++++++++++++++++++++++++++++++
 tests/elevation.test.ts | 31 ++++++++++++++++++++++++
 2 files changed, 84 insertions(+)
 create mode 100644 src/elevation.ts
 create mode 100644 tests/elevation.test.ts

diff --git a/src/elevation.ts b/src/elevation.ts
new file mode 100644
index 0000000..8191928
--- /dev/null
+++ b/src/elevation.ts
@@ -0,0 +1,53 @@
+export type Point = {
+  x: number;
+  y: number;
+};
+
+function calculateDistance(p1: Point, p2: Point): number {
+  return Math.sqrt((p2.x - p1.x) ** 2 + (p2.y - p1.y) ** 2);
+}
+
+function calculateElevationAngle(elevationDifference: number, distance: number): number {
+  return Math.atan2(elevationDifference, distance);
+}
+
+/**
+ * Calculates the maximum heights visible from an observer in a set of directions.
+ * Returns a 2D array. The first dimension is over all Z-values (heights) of the observer, the second dimension
+ * is over the number of theta directions defined with numDirections.
+ * @param grid 2D array of elevations, in the current implementation it needs to be a NxN grid where N is uneven!
+ * @param observer Point of interest, given as the indexes from the grid: [10,10] means that grid[10][10] is the point for
+ * which the elevation angles are calculated.
+ * @param observerZ List of height values of the observer. This will replace the grid[observer] values.
+ * @param numDirections Length of the returned list.
+ * @returns
+ */
+export function getMaxElevationAngles(
+  grid: number[][],
+  observer: Point,
+  observerZ: number[],
+  numDirections: number = 360,
+): number[][] {
+  const maxAngles: number[][] = Array.from({ length: observerZ.length }, () => Array(numDirections).fill(-Infinity));
+  const numRows = grid.length;
+  const numCols = grid[0].length;
+
+  for (let row = 0; row < numRows; row++) {
+    for (let col = 0; col < numCols; col++) {
+      if (row === observer.y && col === observer.x) continue;
+      const targetPoint: Point = { x: col, y: row };
+      const distance = calculateDistance(observer, targetPoint);
+      const angleToTarget = (Math.atan2(targetPoint.y - observer.y, targetPoint.x - observer.x) * 180) / Math.PI;
+      const adjustedAngle = angleToTarget < 0 ? angleToTarget + 360 : angleToTarget;
+      const thetaIndex = Math.floor(adjustedAngle / (360 / numDirections));
+      for (let zIndex = 0; zIndex < observerZ.length; zIndex++) {
+        let elevationDifference = grid[row][col] - observerZ[zIndex];
+        let elevationAngle = calculateElevationAngle(elevationDifference, distance);
+        if (elevationAngle > maxAngles[zIndex][thetaIndex]) {
+          maxAngles[zIndex][thetaIndex] = elevationAngle;
+        }
+      }
+    }
+  }
+  return maxAngles;
+}
diff --git a/tests/elevation.test.ts b/tests/elevation.test.ts
new file mode 100644
index 0000000..cd747ef
--- /dev/null
+++ b/tests/elevation.test.ts
@@ -0,0 +1,31 @@
+import { describe, expect, test } from 'vitest';
+import * as elevation from '../src/elevation';
+
+const observer: elevation.Point = { x: 1, y: 1 }; // Midpoint of the grid
+describe('Test functionalities from elevation.ts: ', () => {
+  const N = 50;
+  const zValuesObserver = [0, 1, 5];
+  const grid = [
+    [1, 0, 2],
+    [3, 0, 4],
+    [5, 0, 6],
+  ];
+
+  const elevationAnglesExpectation = [
+    [1.326, 1.339, 0, 1.295, 1.249, 0.615, 0, 0.955],
+    [1.249, 1.295, -0.785, 1.231, 1.107, 0, -0.785, 0.615],
+    [-0.785, 0.615, -1.373, 0, -1.107, -1.231, -1.373, -1.13],
+  ];
+  const elevationAngles = elevation.getMaxElevationAngles(grid, observer, zValuesObserver, 8);
+  test('Get a list.', () => {
+    const tolerance = 0.01;
+    const allClose =
+      elevationAngles.length === elevationAnglesExpectation.length &&
+      elevationAngles.every(
+        (row, rowIndex) =>
+          row.length === elevationAnglesExpectation[rowIndex].length &&
+          row.every((value, colIndex) => Math.abs(value - elevationAnglesExpectation[rowIndex][colIndex]) <= tolerance),
+      );
+    expect(allClose).toBe(true);
+  });
+});

From f9f2ecad5997eca79c20c61c93eb3a77718a54a8 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 22 Apr 2024 14:21:05 +0200
Subject: [PATCH 03/10] Receive elevation raster in shading simulation #5

---
 src/index.ts | 24 ++++++++++++++++++++----
 1 file changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/index.ts b/src/index.ts
index 5118fdf..6c56e02 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -2,9 +2,9 @@ import * as THREE from 'three';
 import { BufferAttribute, BufferGeometry, TypedArray } from 'three';
 import * as BufferGeometryUtils from 'three/examples/jsm/utils/BufferGeometryUtils.js';
 import { viridis } from './colormaps';
+import * as elevation from './elevation';
 import { getRandomSunVectors } from './sun';
 import * as triangleUtils from './triangleUtils.js';
-import { ArrayType, Triangle } from './triangleUtils.js';
 
 // @ts-ignore
 import { rayTracingWebGL } from './rayTracingWebGL.js';
@@ -20,7 +20,8 @@ import { rayTracingWebGL } from './rayTracingWebGL.js';
 export default class Scene {
   simulationGeometries: Array<BufferGeometry>;
   shadingGeometries: Array<BufferGeometry>;
-  elevationRasters: Array<Array<number[]>>;
+  elevationRasters: Array<number[]>;
+  elevationRasterMidpoint: elevation.Point;
   latitude: number;
   longitude: number;
 
@@ -33,6 +34,7 @@ export default class Scene {
     this.simulationGeometries = [];
     this.shadingGeometries = [];
     this.elevationRasters = [];
+    this.elevationRasterMidpoint = { x: 0, y: 0 };
     this.latitude = latitude;
     this.longitude = longitude;
   }
@@ -59,8 +61,9 @@ export default class Scene {
     this.shadingGeometries.push(geometry);
   }
 
-  addElevationRaster(raster: number[][]) {
-    this.elevationRasters.push(raster);
+  addElevationRaster(raster: number[][], midpoint: elevation.Point) {
+    this.elevationRasters = raster;
+    this.elevationRasterMidpoint = midpoint;
   }
 
   /** @ignore */
@@ -104,6 +107,19 @@ export default class Scene {
     console.log('Simulation package was called to calculate');
     let simulationGeometry = BufferGeometryUtils.mergeGeometries(this.simulationGeometries);
     let shadingGeometry = BufferGeometryUtils.mergeGeometries(this.shadingGeometries);
+    if (this.elevationRasters.length > 0) {
+      let shadingElevationAngles = elevation.getMaxElevationAngles(
+        this.elevationRasters,
+        {
+          x: Math.round(this.elevationRasters.length - 1 / 2),
+          y: Math.round(this.elevationRasters[0].length - 1 / 2),
+        },
+        [this.elevationRasters[this.elevationRasterMidpoint.x][this.elevationRasterMidpoint.y]],
+        // TODO: Right now this only includes one z value taken from the Midpoint, this can be extended
+        // to multiple z points
+      );
+    }
+
     // TODO: This breaks everything, why?
     simulationGeometry = this.refineMesh(simulationGeometry, 0.5); // TODO: make configurable
 

From 854281e03caafde14391f14f9c88e6769ebb6e8f Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 22 Apr 2024 15:23:00 +0200
Subject: [PATCH 04/10] Include shading from elevation to raytracing function
 #5

---
 src/index.ts           | 37 +++++++++++++++++++++++++------------
 src/rayTracingWebGL.ts | 15 +++++++++++----
 src/sun.ts             | 32 +++++++++++++++++++++++++-------
 3 files changed, 61 insertions(+), 23 deletions(-)

diff --git a/src/index.ts b/src/index.ts
index 6c56e02..4361d43 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -107,18 +107,6 @@ export default class Scene {
     console.log('Simulation package was called to calculate');
     let simulationGeometry = BufferGeometryUtils.mergeGeometries(this.simulationGeometries);
     let shadingGeometry = BufferGeometryUtils.mergeGeometries(this.shadingGeometries);
-    if (this.elevationRasters.length > 0) {
-      let shadingElevationAngles = elevation.getMaxElevationAngles(
-        this.elevationRasters,
-        {
-          x: Math.round(this.elevationRasters.length - 1 / 2),
-          y: Math.round(this.elevationRasters[0].length - 1 / 2),
-        },
-        [this.elevationRasters[this.elevationRasterMidpoint.x][this.elevationRasterMidpoint.y]],
-        // TODO: Right now this only includes one z value taken from the Midpoint, this can be extended
-        // to multiple z points
-      );
-    }
 
     // TODO: This breaks everything, why?
     simulationGeometry = this.refineMesh(simulationGeometry, 0.5); // TODO: make configurable
@@ -210,6 +198,31 @@ export default class Scene {
    */
   async rayTrace(midpoints: Float32Array, normals: TypedArray, meshArray: Float32Array, numberSimulations: number) {
     let sunDirections = getRandomSunVectors(numberSimulations, this.latitude, this.longitude);
+    if (this.elevationRasters.length > 0) {
+      let shadingElevationAngles = elevation.getMaxElevationAngles(
+        this.elevationRasters,
+        {
+          x: Math.round(this.elevationRasters.length - 1 / 2),
+          y: Math.round(this.elevationRasters[0].length - 1 / 2),
+        },
+        [this.elevationRasters[this.elevationRasterMidpoint.x][this.elevationRasterMidpoint.y]],
+        // TODO: Right now this only includes one z value taken from the Midpoint, this can be extended
+        // to multiple z points
+      );
+      for (let i = 0; i < sunDirections.spherical.length; i += 2) {
+        const shadingElevationAnglesIndex = Math.round(
+          (sunDirections.spherical[i + 1] / 2 / Math.PI) * shadingElevationAngles.length,
+        );
+        if (shadingElevationAngles[0][shadingElevationAnglesIndex] >= sunDirections.spherical[i]) {
+          sunDirections.cartesian = new Float32Array([
+            ...sunDirections.cartesian.slice(0, i),
+            ...sunDirections.cartesian.slice(i + 3),
+          ]);
+          console.log('One ray was shaded by a nearby mountain.');
+        }
+      }
+    }
+
     return rayTracingWebGL(midpoints, normals, meshArray, sunDirections);
   }
 }
diff --git a/src/rayTracingWebGL.ts b/src/rayTracingWebGL.ts
index bdf9798..d2ea5a3 100644
--- a/src/rayTracingWebGL.ts
+++ b/src/rayTracingWebGL.ts
@@ -10,7 +10,10 @@ export function rayTracingWebGL(
   pointsArray: TypedArray,
   normals: TypedArray,
   trianglesArray: TypedArray,
-  sunDirections: Float32Array,
+  sunDirections: {
+    cartesian: Float32Array;
+    spherical: Float32Array;
+  },
 ): Float32Array | null {
   const N_TRIANGLES = trianglesArray.length / 9;
   const width = pointsArray.length / 3; // Change this to the number of horizontal points in the grid
@@ -178,11 +181,15 @@ export function rayTracingWebGL(
 
   var colorCodedArray = null;
   var isShadowedArray = null;
-  for (var i = 0; i < sunDirections.length; i += 3) {
-    console.log('Simulating sun position #', i / 3, '/', sunDirections.length / 3);
+  for (var i = 0; i < sunDirections.cartesian.length; i += 3) {
+    console.log('Simulating sun position #', i / 3, '/', sunDirections.cartesian.length / 3);
     // TODO: Iterate over sunDirection
     let sunDirectionUniformLocation = gl.getUniformLocation(program, 'u_sun_direction');
-    gl.uniform3fv(sunDirectionUniformLocation, [sunDirections[i], sunDirections[i + 1], sunDirections[i + 2]]);
+    gl.uniform3fv(sunDirectionUniformLocation, [
+      sunDirections.cartesian[i],
+      sunDirections.cartesian[i + 1],
+      sunDirections.cartesian[i + 2],
+    ]);
 
     drawArraysWithTransformFeedback(gl, tf, gl.POINTS, N_POINTS);
 
diff --git a/src/sun.ts b/src/sun.ts
index f6b41cf..71ed7cd 100644
--- a/src/sun.ts
+++ b/src/sun.ts
@@ -1,21 +1,39 @@
 import { getPosition } from 'suncalc';
 
-export function getRandomSunVectors(Ndates: number, lat: number, lon: number): Float32Array {
+/**
+ * Creates arrays of sun vectors. "cartesian" is a vector of length 3*Ndates where every three entries make up one vector.
+ * "spherical" is a vector of length 2*Ndates, where pairs of entries are altitude, azimuth.
+ * @param Ndates
+ * @param lat
+ * @param lon
+ * @returns
+ */
+export function getRandomSunVectors(
+  Ndates: number,
+  lat: number,
+  lon: number,
+): {
+  cartesian: Float32Array;
+  spherical: Float32Array;
+} {
   const sunVectors = new Float32Array(Ndates * 3);
+  const sunVectorsSpherical = new Float32Array(Ndates * 2);
   var i = 0;
   while (i < Ndates) {
     let date = getRandomDate(new Date(2023, 1, 1), new Date(2023, 12, 31));
 
-    const pos = getPosition(date, lat, lon);
-    if (pos.altitude < 0.1 || pos.altitude == Number.NaN) {
+    const posSperical = getPosition(date, lat, lon);
+    if (posSperical.altitude < 0.1 || posSperical.altitude == Number.NaN) {
       continue;
     }
-    sunVectors[3 * i] = -Math.cos(pos.altitude) * Math.sin(pos.azimuth);
-    sunVectors[3 * i + 1] = -Math.cos(pos.altitude) * Math.cos(pos.azimuth);
-    sunVectors[3 * i + 2] = Math.sin(pos.altitude);
+    sunVectors[3 * i] = -Math.cos(posSperical.altitude) * Math.sin(posSperical.azimuth);
+    sunVectors[3 * i + 1] = -Math.cos(posSperical.altitude) * Math.cos(posSperical.azimuth);
+    sunVectors[3 * i + 2] = Math.sin(posSperical.altitude);
+    sunVectorsSpherical[2 * i] = posSperical.altitude;
+    sunVectorsSpherical[2 * i + 1] = posSperical.azimuth;
     i += 1;
   }
-  return sunVectors;
+  return { cartesian: sunVectors, spherical: sunVectorsSpherical };
 }
 
 function getRandomDate(start: Date, end: Date): Date {

From 99785c3a4de32faf19af783744307e5c0495f631 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 22 Apr 2024 15:37:06 +0200
Subject: [PATCH 05/10] Add separation of spherical and cartesian coordinates
 #5

---
 tests/sun.test.ts | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/tests/sun.test.ts b/tests/sun.test.ts
index bcfa57b..cd6c651 100644
--- a/tests/sun.test.ts
+++ b/tests/sun.test.ts
@@ -4,19 +4,25 @@ import * as sun from '../src/sun';
 describe('Test functionalities from sun.ts: ', () => {
   const N = 50;
   let vectors = sun.getRandomSunVectors(N, 0, 0);
-  test('Get Correct number of positions.', () => {
-    expect(vectors.length).toStrictEqual(3 * N);
+  test('Get Correct number of positions for cartesian coordiantes.', () => {
+    expect(vectors.cartesian.length).toStrictEqual(3 * N);
+  });
+  test('Get Correct number of positions for spherical coordiantes.', () => {
+    expect(vectors.spherical.length).toStrictEqual(2 * N);
   });
   test('Get normalized sun vectors.', () => {
     for (let i = 0; i < N / 3; i++) {
-      let length = vectors[3 * i] ** 2 + vectors[3 * i + 1] ** 2 + vectors[3 * i + 2] ** 2;
+      let length = vectors.cartesian[3 * i] ** 2 + vectors.cartesian[3 * i + 1] ** 2 + vectors.cartesian[3 * i + 2] ** 2;
       expect(length).to.closeTo(1, 0.001);
     }
   });
   test('Sun is always above the horizon.', () => {
     for (let i = 0; i < N / 3; i++) {
-      let z = vectors[3 * i + 2];
+      let z = vectors.cartesian[3 * i + 2];
+      let altitude = vectors.spherical[2 * i];
       expect(z).toBeGreaterThan(0);
+      expect(altitude).toBeGreaterThan(0);
+      expect;
     }
   });
 });

From 554d80957aba20ab5fa585eef4d2c626b805ff57 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 22 Apr 2024 17:36:21 +0200
Subject: [PATCH 06/10] Delete unused variable #5

---
 tests/sun.test.ts | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/sun.test.ts b/tests/sun.test.ts
index cd6c651..03bf7d7 100644
--- a/tests/sun.test.ts
+++ b/tests/sun.test.ts
@@ -22,7 +22,6 @@ describe('Test functionalities from sun.ts: ', () => {
       let altitude = vectors.spherical[2 * i];
       expect(z).toBeGreaterThan(0);
       expect(altitude).toBeGreaterThan(0);
-      expect;
     }
   });
 });

From d6e3d858e0148c1fd69e7e981f1c2717e42bad68 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Mon, 6 May 2024 17:54:03 +0200
Subject: [PATCH 07/10] Change 2d array to list of sphericalPoints #5

By doing so, the strange dependency on grid indeces got removed. Now
the list of points has the information of precise locations.
---
 src/elevation.ts        | 101 +++++++++++++++++++++++++++-------------
 tests/elevation.test.ts |  80 +++++++++++++++++++++----------
 2 files changed, 123 insertions(+), 58 deletions(-)

diff --git a/src/elevation.ts b/src/elevation.ts
index 8191928..ab248a4 100644
--- a/src/elevation.ts
+++ b/src/elevation.ts
@@ -1,53 +1,88 @@
 export type Point = {
   x: number;
   y: number;
+  z: number;
 };
 
-function calculateDistance(p1: Point, p2: Point): number {
-  return Math.sqrt((p2.x - p1.x) ** 2 + (p2.y - p1.y) ** 2);
+export type SphericalPoint = {
+  altitude: number;
+  azimuth: number;
+};
+
+export function fillMissingAltitudes(maxAngles: SphericalPoint[]): void {
+  // First copy the maxAngles to a newAngles list, so that changes
+  // in the list do not affect the algorithm
+  let newAngles = maxAngles.map((angle) => ({ ...angle }));
+  for (let i = 0; i < newAngles.length; i++) {
+    if (newAngles[i].altitude != -Infinity) {
+      continue;
+    }
+    let distance = 1;
+    while (true) {
+      let prevIndex = (i - distance + newAngles.length) % newAngles.length;
+      let nextIndex = (i + distance) % newAngles.length;
+
+      if (maxAngles[nextIndex].altitude !== -Infinity) {
+        newAngles[i].altitude = maxAngles[nextIndex].altitude;
+        break;
+      } else if (maxAngles[prevIndex].altitude !== -Infinity) {
+        newAngles[i].altitude = maxAngles[prevIndex].altitude;
+        break;
+      } else distance++;
+    }
+  }
+  // Overwrite the maxAngles to make changes in this vector global
+  for (let i = 0; i < maxAngles.length; i++) {
+    maxAngles[i] = newAngles[i];
+  }
 }
 
-function calculateElevationAngle(elevationDifference: number, distance: number): number {
-  return Math.atan2(elevationDifference, distance);
+/**
+ *
+ * @param start
+ * @param end
+ * @returns azimuth from 0 to 2*PI and altitude from 0 to PI/2, where altitude = 0 is facing directly upwards
+ */
+export function calculateSphericalCoordinates(start: Point, end: Point): { azimuth: number; altitude: number } {
+  const dx = end.x - start.x;
+  const dy = end.y - start.y;
+  const dz = end.z - start.z;
+
+  const r = Math.sqrt(dx * dx + dy * dy + dz * dz);
+  const altitude = Math.acos(dz / r);
+  let azimuth = Math.atan2(dy, dx);
+
+  if (azimuth < 0) {
+    azimuth += 2 * Math.PI; // Adjust azimuth to be from 0 to 2PI
+  }
+
+  return { azimuth, altitude };
 }
 
 /**
  * Calculates the maximum heights visible from an observer in a set of directions.
- * Returns a 2D array. The first dimension is over all Z-values (heights) of the observer, the second dimension
- * is over the number of theta directions defined with numDirections.
- * @param grid 2D array of elevations, in the current implementation it needs to be a NxN grid where N is uneven!
+ * Returns a list of spherical points of length numDirections.
+ * @param elevation list of points with x,y,z component
  * @param observer Point of interest, given as the indexes from the grid: [10,10] means that grid[10][10] is the point for
  * which the elevation angles are calculated.
- * @param observerZ List of height values of the observer. This will replace the grid[observer] values.
  * @param numDirections Length of the returned list.
  * @returns
  */
-export function getMaxElevationAngles(
-  grid: number[][],
-  observer: Point,
-  observerZ: number[],
-  numDirections: number = 360,
-): number[][] {
-  const maxAngles: number[][] = Array.from({ length: observerZ.length }, () => Array(numDirections).fill(-Infinity));
-  const numRows = grid.length;
-  const numCols = grid[0].length;
-
-  for (let row = 0; row < numRows; row++) {
-    for (let col = 0; col < numCols; col++) {
-      if (row === observer.y && col === observer.x) continue;
-      const targetPoint: Point = { x: col, y: row };
-      const distance = calculateDistance(observer, targetPoint);
-      const angleToTarget = (Math.atan2(targetPoint.y - observer.y, targetPoint.x - observer.x) * 180) / Math.PI;
-      const adjustedAngle = angleToTarget < 0 ? angleToTarget + 360 : angleToTarget;
-      const thetaIndex = Math.floor(adjustedAngle / (360 / numDirections));
-      for (let zIndex = 0; zIndex < observerZ.length; zIndex++) {
-        let elevationDifference = grid[row][col] - observerZ[zIndex];
-        let elevationAngle = calculateElevationAngle(elevationDifference, distance);
-        if (elevationAngle > maxAngles[zIndex][thetaIndex]) {
-          maxAngles[zIndex][thetaIndex] = elevationAngle;
-        }
-      }
+export function getMaxElevationAngles(elevation: Point[], observer: Point, numDirections: number = 360): SphericalPoint[] {
+  let maxAngles: SphericalPoint[] = Array.from({ length: numDirections }, (_, index) => ({
+    azimuth: index * ((2 * Math.PI) / numDirections),
+    altitude: -Infinity,
+  }));
+
+  for (let point of elevation) {
+    const { azimuth, altitude } = calculateSphericalCoordinates(observer, point);
+    console.log(azimuth, altitude);
+    const closestIndex = Math.round(azimuth / ((2 * Math.PI) / numDirections)) % numDirections;
+
+    if (altitude > maxAngles[closestIndex].altitude) {
+      maxAngles[closestIndex].altitude = altitude;
     }
   }
+  fillMissingAltitudes(maxAngles);
   return maxAngles;
 }
diff --git a/tests/elevation.test.ts b/tests/elevation.test.ts
index cd747ef..6bb0715 100644
--- a/tests/elevation.test.ts
+++ b/tests/elevation.test.ts
@@ -1,31 +1,61 @@
 import { describe, expect, test } from 'vitest';
 import * as elevation from '../src/elevation';
 
-const observer: elevation.Point = { x: 1, y: 1 }; // Midpoint of the grid
-describe('Test functionalities from elevation.ts: ', () => {
-  const N = 50;
-  const zValuesObserver = [0, 1, 5];
-  const grid = [
-    [1, 0, 2],
-    [3, 0, 4],
-    [5, 0, 6],
-  ];
+describe('calculateSphericalCoordinates', () => {
+  test('should calculate the correct spherical coordinates', () => {
+    const start = { x: 0, y: 0, z: 0 };
+    const ends = [
+      { x: 0, y: 0, z: 1 },
+      { x: 1, y: 0, z: 1 },
+      { x: 0, y: -1, z: 1 },
+    ];
+    const expectedResults = [
+      { altitude: 0, azimuth: 0 },
+      { altitude: Math.PI / 4, azimuth: 0 },
+      { altitude: Math.PI / 4, azimuth: (3 / 2) * Math.PI },
+    ];
 
-  const elevationAnglesExpectation = [
-    [1.326, 1.339, 0, 1.295, 1.249, 0.615, 0, 0.955],
-    [1.249, 1.295, -0.785, 1.231, 1.107, 0, -0.785, 0.615],
-    [-0.785, 0.615, -1.373, 0, -1.107, -1.231, -1.373, -1.13],
-  ];
-  const elevationAngles = elevation.getMaxElevationAngles(grid, observer, zValuesObserver, 8);
-  test('Get a list.', () => {
-    const tolerance = 0.01;
-    const allClose =
-      elevationAngles.length === elevationAnglesExpectation.length &&
-      elevationAngles.every(
-        (row, rowIndex) =>
-          row.length === elevationAnglesExpectation[rowIndex].length &&
-          row.every((value, colIndex) => Math.abs(value - elevationAnglesExpectation[rowIndex][colIndex]) <= tolerance),
-      );
-    expect(allClose).toBe(true);
+    ends.forEach((end, index) => {
+      const result = elevation.calculateSphericalCoordinates(start, end);
+      const expected = expectedResults[index];
+      expect(result.azimuth).toBeCloseTo(expected.azimuth);
+      expect(result.altitude).toBeCloseTo(expected.altitude);
+    });
+  });
+});
+
+describe('fillMissingAltitudes', () => {
+  test('should fill negative infinity altitude with the nearest non-negative infinity altitude', () => {
+    const points: elevation.SphericalPoint[] = [
+      { altitude: -Infinity, azimuth: 0 },
+      { altitude: 10, azimuth: 90 },
+      { altitude: -Infinity, azimuth: 180 },
+      { altitude: -Infinity, azimuth: 230 },
+      { altitude: -Infinity, azimuth: 240 },
+      { altitude: 20, azimuth: 270 },
+    ];
+
+    elevation.fillMissingAltitudes(points);
+
+    expect(points[0].altitude).toBe(10);
+    expect(points[2].altitude).toBe(10);
+    expect(points[3].altitude).toBe(20);
+    expect(points[4].altitude).toBe(20);
+  });
+});
+
+describe('getMaxElevationAngles', () => {
+  test('should correctly calculate the maximum elevation angles for given elevation points and observer', () => {
+    const elevations: elevation.Point[] = [
+      { x: 1, y: 1, z: 2 },
+      { x: 1, y: -1, z: 4 },
+      { x: -1, y: -1, z: 6 },
+      { x: -1, y: 1, z: 8 },
+    ];
+    const observer: elevation.Point = { x: 0, y: 0, z: 0 };
+    const numDirections = 20;
+    const result: elevation.SphericalPoint[] = elevation.getMaxElevationAngles(elevations, observer, numDirections);
+    console.log(result);
+    expect(result).to.be.an('array').that.has.lengthOf(numDirections);
   });
 });

From 3ac9980b817bf6ba0f2305bceace7e84e49eb3ec Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Tue, 7 May 2024 08:01:29 +0200
Subject: [PATCH 08/10] Change type of elevation to list of points #5

---
 src/index.ts | 37 ++++++++++---------------------------
 1 file changed, 10 insertions(+), 27 deletions(-)

diff --git a/src/index.ts b/src/index.ts
index 4361d43..6a11b6c 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -20,7 +20,7 @@ import { rayTracingWebGL } from './rayTracingWebGL.js';
 export default class Scene {
   simulationGeometries: Array<BufferGeometry>;
   shadingGeometries: Array<BufferGeometry>;
-  elevationRasters: Array<number[]>;
+  elevationRaster: Array<elevation.Point>;
   elevationRasterMidpoint: elevation.Point;
   latitude: number;
   longitude: number;
@@ -33,8 +33,8 @@ export default class Scene {
   constructor(latitude: number, longitude: number) {
     this.simulationGeometries = [];
     this.shadingGeometries = [];
-    this.elevationRasters = [];
-    this.elevationRasterMidpoint = { x: 0, y: 0 };
+    this.elevationRaster = [];
+    this.elevationRasterMidpoint = { x: 0, y: 0, z: 0 };
     this.latitude = latitude;
     this.longitude = longitude;
   }
@@ -61,8 +61,8 @@ export default class Scene {
     this.shadingGeometries.push(geometry);
   }
 
-  addElevationRaster(raster: number[][], midpoint: elevation.Point) {
-    this.elevationRasters = raster;
+  addElevationRaster(raster: elevation.Point[], midpoint: elevation.Point) {
+    this.elevationRaster = raster;
     this.elevationRasterMidpoint = midpoint;
   }
 
@@ -198,31 +198,14 @@ export default class Scene {
    */
   async rayTrace(midpoints: Float32Array, normals: TypedArray, meshArray: Float32Array, numberSimulations: number) {
     let sunDirections = getRandomSunVectors(numberSimulations, this.latitude, this.longitude);
-    if (this.elevationRasters.length > 0) {
+    if (this.elevationRaster.length > 0) {
       let shadingElevationAngles = elevation.getMaxElevationAngles(
-        this.elevationRasters,
-        {
-          x: Math.round(this.elevationRasters.length - 1 / 2),
-          y: Math.round(this.elevationRasters[0].length - 1 / 2),
-        },
-        [this.elevationRasters[this.elevationRasterMidpoint.x][this.elevationRasterMidpoint.y]],
-        // TODO: Right now this only includes one z value taken from the Midpoint, this can be extended
-        // to multiple z points
+        this.elevationRaster,
+        this.elevationRasterMidpoint,
+        numberSimulations,
       );
-      for (let i = 0; i < sunDirections.spherical.length; i += 2) {
-        const shadingElevationAnglesIndex = Math.round(
-          (sunDirections.spherical[i + 1] / 2 / Math.PI) * shadingElevationAngles.length,
-        );
-        if (shadingElevationAngles[0][shadingElevationAnglesIndex] >= sunDirections.spherical[i]) {
-          sunDirections.cartesian = new Float32Array([
-            ...sunDirections.cartesian.slice(0, i),
-            ...sunDirections.cartesian.slice(i + 3),
-          ]);
-          console.log('One ray was shaded by a nearby mountain.');
-        }
-      }
     }
-
+    //TODO: add shading of elevation here
     return rayTracingWebGL(midpoints, normals, meshArray, sunDirections);
   }
 }

From 67e84883c33b7e3ad578e1b1d2e3506e4eb0a795 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Tue, 7 May 2024 08:10:01 +0200
Subject: [PATCH 09/10] Correct numDirections parameter #5

---
 src/index.ts | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/index.ts b/src/index.ts
index 6a11b6c..0eda68e 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -193,6 +193,7 @@ export default class Scene {
    * @param midpoints midpoints of triangles for which to calculate intensities
    * @param normals normals for each midpoint
    * @param meshArray array of vertices for the shading mesh
+   * @param numberSimulations
    * @return
    * @memberof Scene
    */
@@ -202,7 +203,7 @@ export default class Scene {
       let shadingElevationAngles = elevation.getMaxElevationAngles(
         this.elevationRaster,
         this.elevationRasterMidpoint,
-        numberSimulations,
+        sunDirections.spherical.length / 2,
       );
     }
     //TODO: add shading of elevation here

From 6f1a9a4c2d87250150881c575061148cc1285a03 Mon Sep 17 00:00:00 2001
From: Florian Kotthoff <74312290+FlorianK13@users.noreply.github.com>
Date: Tue, 7 May 2024 08:10:16 +0200
Subject: [PATCH 10/10] Correct function description #5

---
 src/elevation.ts | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/elevation.ts b/src/elevation.ts
index ab248a4..7322f8a 100644
--- a/src/elevation.ts
+++ b/src/elevation.ts
@@ -63,9 +63,8 @@ export function calculateSphericalCoordinates(start: Point, end: Point): { azimu
  * Calculates the maximum heights visible from an observer in a set of directions.
  * Returns a list of spherical points of length numDirections.
  * @param elevation list of points with x,y,z component
- * @param observer Point of interest, given as the indexes from the grid: [10,10] means that grid[10][10] is the point for
- * which the elevation angles are calculated.
- * @param numDirections Length of the returned list.
+ * @param observer Point of interest for which the elevation angles are calculated.
+ * @param numDirections Number of steps for the azimuth angle.
  * @returns
  */
 export function getMaxElevationAngles(elevation: Point[], observer: Point, numDirections: number = 360): SphericalPoint[] {