diff --git a/CHANGELOG.md b/CHANGELOG.md
index e4f97fa..ea3e031 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to
 [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## 4.0.0 - 2021-08-25
+
+### Added
+
+- Option to show 22r and 46r guide circles, horizon line, nadir and zenith
+
+### Changed
+
+- Crystal geometry generation logic was rewritten to keep pyramid apex angle
+  constant regardless of varying prism face distances to the crystal C-axis
+
 ## 3.3.0 - 2021-05-07
 
 ### Changed
diff --git a/README.md b/README.md
index 04fba8b..b4f690d 100644
--- a/README.md
+++ b/README.md
@@ -138,6 +138,7 @@ These settings affect how the results of the simulation are shown on the screen.
 - **Brightness:** Alters the total brightness of the image, much like an exposure adjustment on cameras
 - **Hide sub-horizon:** Hides any halos below the horizon level
 - **Lock to light source:** Locks the camera to the sun
+- **Show guides:** Draws markings for horizon, zenith, nadir, 22r and 46r circles
 
 ### Atmosphere settings
 
@@ -238,4 +239,5 @@ problems like this are encountered.
 - [Lauri Kangas](https://github.com/lkangas) for providing tons of reading material and debugging help
 - [Panu Lahtinen](https://github.com/pnuu) for additional Linux support
 - Jukka Ruoskanen for making HaloPoint 2.0 back in the day and inspiring me to start working on HaloRay
+- Marko Riikonen for giving valuable feedback and reporting bugs
 - [Jaakko Lehtinen](https://users.aalto.fi/~lehtinj7/) for super valuable lessons in computer graphics
diff --git a/appveyor.yml b/appveyor.yml
index 44e1799..5711e17 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -1,4 +1,4 @@
-version: "3.3.0-{build}"
+version: "4.0.0-{build}"
 branches:
   only:
     - master
diff --git a/src/haloray-core/gui/crystalPreview/previewRenderArea.cpp b/src/haloray-core/gui/crystalPreview/previewRenderArea.cpp
index e1ddabb..80a8bc9 100644
--- a/src/haloray-core/gui/crystalPreview/previewRenderArea.cpp
+++ b/src/haloray-core/gui/crystalPreview/previewRenderArea.cpp
@@ -39,44 +39,46 @@ QSize PreviewRenderArea::sizeHint() const
 void PreviewRenderArea::paintEvent(QPaintEvent *)
 {
     const int numVertices = 24;
-    initializeGeometry(m_vertices, numVertices);
-    float largestDimension = getFurthestVertexDistance(m_vertices, numVertices);
+    initializeGeometry(m_vertices);
 
-    QMatrix4x4 transformMat;
-    transformMat.scale(600.0f);
-    transformMat.perspective(90.0f, 1.0, 0.01f, 5.0f);
-    transformMat.lookAt(QVector3D(-2.0f, 2.0f, 2.0f), QVector3D(0.0f, 0.0f, 0.0f), QVector3D(0.0f, 1.0f, 0.0f));
-    transformMat.scale(1.0f / largestDimension);
+    QMatrix4x4 perspectiveMat;
+    perspectiveMat.perspective(90.0f, 1.0, 0.01f, 100.0f);
 
-    QPainter painter(this);
-    painter.setRenderHint(QPainter::RenderHint::Antialiasing);
-    QPen pen;
-    pen.setColor(QColor(0, 0, 0));
-    pen.setWidth(3);
-    painter.setPen(pen);
-
-    painter.rotate(180.0);
-    painter.translate(-width() / 2, -height() / 2);
-    int side = qMin(width(), height());
-    painter.scale(side / 500.0f, side / 500.0f);
+    QMatrix4x4 viewMat;
+    viewMat.lookAt(QVector3D(5.0f, 5.0f, 5.0f), QVector3D(0.0f, 0.0f, 0.0f), QVector3D(0.0f, 1.0f, 0.0f));
 
-    QMatrix4x4 orientationMatrix = getCrystalOrientationMatrix();
+    float largestDimension = getFurthestVertexDistance(m_vertices, numVertices);
+    QMatrix4x4 orientationMat = getCrystalOrientationMatrix();
+    QMatrix4x4 modelMat;
+    modelMat.scale(1.0f / largestDimension);
 
     QVector4D mappedVertices[numVertices];
     for (int i = 0; i < numVertices; ++i)
     {
-        mappedVertices[i] = transformMat * orientationMatrix * QVector4D(m_vertices[i], 1.0f);
+        mappedVertices[i] = perspectiveMat * viewMat * orientationMat * modelMat * QVector4D(m_vertices[i], 1.0f);
     }
 
-    QPoint points[numVertices];
+    QPointF points[numVertices];
     std::transform(
                 mappedVertices,
                 mappedVertices + numVertices,
                 points,
                 [](QVector4D vertex) {
-        return (vertex / vertex.w()).toPoint();
+        return (vertex / vertex.w()).toPointF();
     });
 
+    QPainter painter(this);
+    painter.setRenderHint(QPainter::RenderHint::Antialiasing);
+    painter.translate(width() / 2, height() / 2);
+    int side = qMin(width(), height());
+    painter.scale(side, -side);
+    painter.scale(3.5, 3.5);
+
+    QPen pen;
+    pen.setColor(QColor(0, 0, 0));
+    pen.setWidthF(0.5/side);
+    painter.setPen(pen);
+
     painter.drawPolygon(points, 6);
     painter.drawPolygon(points + 6 , 6);
     painter.drawPolygon(points + 12, 6);
@@ -84,28 +86,126 @@ void PreviewRenderArea::paintEvent(QPaintEvent *)
 
     for (int i = 0; i < 6; ++i)
     {
-        QPolygon edgePoints;
+        QPolygonF edgePoints;
         edgePoints << points[i] << points [i + 6] << points[i + 12] << points[i + 18];
         painter.drawPolyline(edgePoints);
     }
+
+    // Drawing the axis lines is helpful for debugging
+    //drawAxisLines(perspectiveMat * viewMat, &painter);
 }
 
-QVector2D lineIntersect(QVector2D line1, QVector2D line2)
+void PreviewRenderArea::drawAxisLines(QMatrix4x4 viewMat, QPainter *painter)
 {
-    float p1 = line1.x();
-    float theta1 = line1.y();
+    auto axisLength = 10.0f;
+    auto origin = viewMat * QVector4D(0.0f, 0.0f, 0.0f, 1.0f);
 
-    float p2 = line2.x();
-    float theta2 = line2.y();
+    QPen pen;
+    pen.setWidth(0);
+    pen.setColor(QColor(255, 0, 0));
+    painter->setPen(pen);
+    auto xAxis = viewMat * QVector4D(axisLength, 0.0f, 0.0f, 1.0f);
+    painter->drawLine((origin / origin.w()).toPoint(), (xAxis / xAxis.w()).toPointF());
+
+    pen.setColor(QColor(0, 255, 0));
+    painter->setPen(pen);
+    auto yAxis = viewMat * QVector4D(0.0f, axisLength, 0.0f, 1.0f);
+    painter->drawLine((origin / origin.w()).toPoint(), (yAxis / yAxis.w()).toPointF());
+
+    pen.setColor(QColor(0, 0, 255));
+    painter->setPen(pen);
+    auto zAxis = viewMat * QVector4D(0.0f, 0.0f, axisLength, 1.0f);
+    painter->drawLine((origin / origin.w()).toPoint(), (zAxis / zAxis.w()).toPointF());
+}
+
+int getPrevious(int i)
+{
+    return (((i - 1) % 6) + 6) % 6;
+}
+
+int getNext(int i)
+{
+    return (i + 1) % 6;
+}
+
+QVector2D rotate(double angle, double x, double y)
+{
+    return QVector2D(cos(angle) * x - sin(angle) * y, sin(angle) * x + cos(angle) * y);
+}
+
+float determinant3x3(QMatrix3x3 mat)
+{
+    auto a = mat.constData()[0];
+    auto b = mat.constData()[1];
+    auto c = mat.constData()[2];
+    auto d = mat.constData()[3];
+    auto e = mat.constData()[4];
+    auto f = mat.constData()[5];
+    auto g = mat.constData()[6];
+    auto h = mat.constData()[7];
+    auto i = mat.constData()[8];
+    return a * (e * i - h * f) - b * (d * i - g * f) + c * (d * h - g * e);
+}
+
+void generateApexNormals(double apexAngle, QVector3D *apexNormals)
+{
+    auto halfApex = apexAngle / 2.0;
+    QMatrix4x4 normalTiltMat;
+    normalTiltMat.rotate(-halfApex * 180.0f / PI, 1.0f, 0.0f, 0.0f);
+    for (auto i = 0; i < 6; ++i)
+    {
+        auto rotAngle = i * PI / 3.0;
+        /* Initial normal vector is first tilted around the X axis
+         * according to the apex angle and then rotated around the Y
+         * axis in 60 degree increments. */
+        QMatrix4x4 normalRotationMat;
+        normalRotationMat.rotate(rotAngle * 180.0f / PI, 0.0f, 1.0f, 0.0f);
+        apexNormals[i] = normalRotationMat * normalTiltMat * QVector3D(0.0f, 0.0f, 1.0f);
+    }
+}
 
-    float deltaSine = sin(theta2 - theta1);
-    float x = (p1 * sin(theta2) - p2 * sin(theta1)) / deltaSine;
-    float y = (p2 * cos(theta1) - p1 * cos(theta2)) / deltaSine;
+float getMaximumApexHeight(QVector3D *normals, QVector3D *vertices, float *prismFaceDistances, float apexAngle, int vertexOffset)
+{
+    float maxApexHeight = std::numeric_limits<float>::max();
+    // Intersect three adjacent pyramid faces to find maximum height of pyramid cap based on degenerating faces
+    for (auto face = 0; face < 6; ++face)
+    {
+        auto prevFace = getPrevious(face);
+        auto nextFace = getNext(face);
+
+        auto nPrev = normals[prevFace];
+        auto nCurr = normals[face];
+        auto nNext = normals[nextFace];
+
+        auto prevVert = vertices[prevFace + vertexOffset];
+        auto nextVert = vertices[face + vertexOffset];
+
+        auto numerator = QVector3D::dotProduct(prevVert, nPrev) * QVector3D::crossProduct(nCurr, nNext) +
+                QVector3D::dotProduct(prevVert, nCurr) * QVector3D::crossProduct(nNext, nPrev) +
+                QVector3D::dotProduct(nextVert, nNext) * QVector3D::crossProduct(nPrev, nCurr);
+        float detMatrixContents[] = {
+            nPrev.x(), nCurr.x(), nNext.x(),
+            nPrev.y(), nCurr.y(), nNext.y(),
+            nPrev.z(), nCurr.z(), nNext.z()
+        };
+        auto detMatrix = QMatrix3x3(detMatrixContents);
+        auto denominator = determinant3x3(detMatrix);
+        auto faceHeight = abs((numerator / denominator).y());
+        maxApexHeight = fminf(faceHeight, maxApexHeight);
+    }
 
-    return QVector2D(x, y);
+    // Find maximum height of pyramid cap based on angle and face distances
+    for (auto i = 0; i < 3; ++i)
+    {
+        auto dist = prismFaceDistances[i];
+        auto distOpposite = prismFaceDistances[i + 3];
+        auto h = (dist + distOpposite) / (2.0 * tan(apexAngle / 2.0));
+        maxApexHeight = fminf(h, maxApexHeight);
+    }
+    return maxApexHeight;
 }
 
-void PreviewRenderArea::initializeGeometry(QVector3D *vertices, int numVertices)
+void PreviewRenderArea::initializeGeometry(QVector3D *vertices)
 {
     float caRatioAverage = getFromModel(m_populationIndex, CrystalModel::CaRatioAverage).toFloat();
     float upperApexAngle = degToRad(getFromModel(m_populationIndex, CrystalModel::UpperApexAngle).toFloat());
@@ -121,87 +221,112 @@ void PreviewRenderArea::initializeGeometry(QVector3D *vertices, int numVertices)
         getFromModel(m_populationIndex, CrystalModel::PrismFaceDistance6).toFloat(),
     };
 
-    float deltaAngle = degToRad(60.0f);
     QVector2D hexagonCorners[6];
-    /* The sqrt(3)/2 multiplier makes the default crystal such
-       that the distance of a vertex from the C axis is 1.0 */
-    float sizeScaler = sqrt(3.0f) / 2.0f;
-    for (int face = 0; face < 6; ++face)
-    {
-        int previousFace = face == 0 ? 5 : face - 1;
-        int nextFace = face == 5 ? 0 : face + 1;
-
-        float previousAngle = (face + 1) * deltaAngle;
-        float currentAngle = previousAngle + deltaAngle;
-        float nextAngle = previousAngle + 2.0 * deltaAngle;
 
-        float previousDistance = sizeScaler * prismFaceDistances[previousFace];
-        float currentDistance = sizeScaler * prismFaceDistances[face];
-        float nextDistance = sizeScaler * prismFaceDistances[nextFace];
-
-        QVector2D previousLine = QVector2D(previousDistance, previousAngle);
-        QVector2D currentLine = QVector2D(currentDistance, currentAngle);
-        QVector2D nextLine = QVector2D(nextDistance, nextAngle);
-
-        QVector2D previousCurrentIntersection = lineIntersect(previousLine, currentLine);
-        QVector2D currentNextIntersection = lineIntersect(currentLine, nextLine);
-        QVector2D previousNextIntersection = lineIntersect(previousLine, nextLine);
-
-        float previousCurrentIntersectionDistance = previousCurrentIntersection.length();
-        float currentNextIntersectionDistance = currentNextIntersection.length();
-        float previousNextIntersectionDistance = previousNextIntersection.length();
-
-        QVector2D v1 = previousCurrentIntersectionDistance < previousNextIntersectionDistance ? previousCurrentIntersection : previousNextIntersection;
-        QVector2D v2 = currentNextIntersectionDistance < previousNextIntersectionDistance ? currentNextIntersection : previousNextIntersection;
+    // Calculate initial hexagon corners
+    for (auto i = 0; i < 6; ++i)
+    {
+        auto d1 = prismFaceDistances[i];
+        auto d2 = prismFaceDistances[getNext(i)];
 
-        if (face > 0 && previousNextIntersectionDistance > hexagonCorners[face].length())
-        {
-            v1 = hexagonCorners[face];
-        }
+        auto angle = -i * PI / 3.0;
+        auto x_stat = 2.0 * d2 / sqrt(3.0) - d1 / sqrt(3.0);
+        auto y_stat = d1;
+        hexagonCorners[i] = rotate(angle, x_stat, y_stat);
+    }
 
-        if (face == 5 && previousNextIntersectionDistance > hexagonCorners[nextFace].length())
+    // Fix denegerate prism faces
+    for (auto face = 0; face < 6; ++face)
+    {
+        auto prevFace = getPrevious(face);
+        auto nextFace = getNext(face);
+        auto angle = -face * PI / 3.0;
+
+        auto d1 = prismFaceDistances[face];
+        auto d2 = prismFaceDistances[nextFace];
+        auto d3 = prismFaceDistances[prevFace];
+        if (d1 > d2 + d3)
         {
-            v2 = hexagonCorners[nextFace];
+            auto x_stat = d2 / sqrt(3.0) - d3 / sqrt(3.0);
+            auto y_stat = d2 + d3;
+            auto rotatedPoint = rotate(angle, x_stat, y_stat);
+            hexagonCorners[face] = rotatedPoint;
+            hexagonCorners[prevFace] = rotatedPoint;
         }
+    }
 
-        hexagonCorners[face] = v1;
-        hexagonCorners[nextFace] = v2;
+    /* Scaling value makes sure eventual A axis length is 2.0, so C/A ratio
+     * can be easily corrected. */
+    auto hexagonScaler = (hexagonCorners[1] - hexagonCorners[4]).length();
+    for (auto i = 0; i < 6; ++i)
+    {
+        hexagonCorners[i] *= 2.0f / hexagonScaler;
     }
 
-    for (int face = 0; face < 6; ++face)
+    for (auto face = 0; face < 6; ++face)
     {
         QVector2D *vertex = &hexagonCorners[face];
         // Corresponding vertices of each crystal layer have the same X and Z coordinates
-        vertices[face] = QVector3D(vertex->x(), 1.0f, vertex->y());
-        vertices[face + 6] = QVector3D(vertex->x(), 1.0f, vertex->y());
-        vertices[face + 12] = QVector3D(vertex->x(), -1.0f, vertex->y());
-        vertices[face + 18] = QVector3D(vertex->x(), -1.0f, vertex->y());
+        // First six vertices are the top apex cap
+        // Next six vertices are the top of the base hexagonal crystal
+        // Next six vertices are the bottom of the base hexagonal crystal
+        // Last six vertices are the bottom apex cap
+        vertices[face] = QVector3D(vertex->x(), 0.0f, vertex->y());
+        vertices[face + 6] = QVector3D(vertex->x(), 0.0f, vertex->y());
+        vertices[face + 12] = QVector3D(vertex->x(), 0.0f, vertex->y());
+        vertices[face + 18] = QVector3D(vertex->x(), 0.0f, vertex->y());
     }
 
-    // Stretch the crystal to correct C/A ratio
-    float caMultiplier = caRatioAverage;
-    for (int i = 0; i < numVertices; ++i)
+    if (upperApexHeightAverage > 0.0 && upperApexAngle < PI && upperApexAngle > 0.0)
     {
-        vertices[i].setY(vertices[i].y() * caMultiplier);
+        // Generate normals for upper pyramid cap
+        QVector3D upperApexNormals[6];
+        generateApexNormals(upperApexAngle, upperApexNormals);
+
+        // Set upper pyramid cap vertex positions
+        float maxUpperApexHeight = getMaximumApexHeight(upperApexNormals, vertices, prismFaceDistances, upperApexAngle, 0);
+        for (auto i = 0; i < 6; ++i)
+        {
+            auto next = getNext(i);
+            auto pyramidEdge = QVector3D::crossProduct(upperApexNormals[i], upperApexNormals[next]);
+            vertices[i] += upperApexHeightAverage * maxUpperApexHeight * pyramidEdge / pyramidEdge.y();
+        }
     }
 
-    // Scale pyramid caps
-    float upperApexMaxHeight = sizeScaler / tan(upperApexAngle / 2.0);
-    float lowerApexMaxHeight = sizeScaler / tan(lowerApexAngle / 2.0);
+    if (lowerApexHeightAverage > 0.0 && lowerApexAngle < PI && lowerApexAngle > 0.0)
+    {
+        // Generate normals for lower pyramid cap
+        QVector3D lowerApexNormals[6];
+        generateApexNormals(lowerApexAngle, lowerApexNormals);
+        for (auto i = 0; i < 6; ++i)
+        {
+            lowerApexNormals[i].setY(-lowerApexNormals[i].y());
+        }
 
-    float upperApexHeight = upperApexHeightAverage;
-    float lowerApexHeight = lowerApexHeightAverage;
+        // Set lower pyramid cap vertex positions
+        float maxLowerApexHeight = getMaximumApexHeight(lowerApexNormals, vertices, prismFaceDistances, lowerApexAngle, 18);
+        for (auto i = 0; i < 6; ++i)
+        {
+            auto next = getNext(i);
+            auto pyramidEdge = QVector3D::crossProduct(lowerApexNormals[i], lowerApexNormals[next]);
+            vertices[i + 18] -= lowerApexHeightAverage * maxLowerApexHeight * pyramidEdge / pyramidEdge.y();
+        }
+    }
 
-    for (int i = 0; i < 6; ++i)
+    // Scale crystal vertically to have correct C/A ratio
+    for (auto i = 0; i < 12; ++i)
     {
-        vertices[i] = QVector3D(
-                    vertices[i].x() * (1.0 - upperApexHeight),
-                    vertices[i].y() + upperApexHeight * upperApexMaxHeight,
-                    vertices[i].z() * (1.0 - upperApexHeight));
-        vertices[numVertices - i - 1] = QVector3D(
-                    vertices[numVertices - i - 1].x() * (1.0 - lowerApexHeight),
-                    vertices[numVertices - i - 1].y() - lowerApexHeight * lowerApexMaxHeight,
-                    vertices[numVertices - i - 1].z() * (1.0 - lowerApexHeight));
+        vertices[i].setY(vertices[i].y() + caRatioAverage);
+        vertices[i + 12].setY(vertices[i + 12].y() - caRatioAverage);
+    }
+
+    // Rotate crystal around C-axis so that face numbering follows conventions
+    // Prism face 0 (Face 3 in the UI) should be up in a column Parry position
+    QMatrix4x4 conventionMatrix;
+    conventionMatrix.rotate(-90.0f, 0.0f, 1.0f, 0.0f);
+    for (auto i = 0; i < 24; ++i)
+    {
+        vertices[i] = conventionMatrix * vertices[i];
     }
 }
 
@@ -221,8 +346,10 @@ QMatrix4x4 PreviewRenderArea::getCrystalOrientationMatrix() const
     float tilt = getFromModel(m_populationIndex, CrystalModel::TiltAverage).toFloat();
     float rotation = getFromModel(m_populationIndex, CrystalModel::RotationAverage).toFloat();
     QMatrix4x4 orientationMatrix;
-    orientationMatrix.rotate(tilt, QVector3D(0.0f, 0.0f, 1.0f));
+    // First rotate around Y and tent tile around Z
+    orientationMatrix.rotate(-tilt, QVector3D(0.0f, 0.0f, 1.0f));
     orientationMatrix.rotate(rotation, QVector3D(0.0f, 1.0f, 0.0f));
+
     return orientationMatrix;
 }
 
diff --git a/src/haloray-core/gui/crystalPreview/previewRenderArea.h b/src/haloray-core/gui/crystalPreview/previewRenderArea.h
index 9e237ba..a9aa172 100644
--- a/src/haloray-core/gui/crystalPreview/previewRenderArea.h
+++ b/src/haloray-core/gui/crystalPreview/previewRenderArea.h
@@ -15,19 +15,20 @@ class PreviewRenderArea : public QWidget
     Q_OBJECT
 public:
     PreviewRenderArea(CrystalModel *crystals, QWidget *parent = nullptr);
+    QSize sizeHint() const override;
 
 public slots:
     void onPopulationSelectionChange(int index);
-    QSize sizeHint() const override;
 
 protected:
     void paintEvent(QPaintEvent *event) override;
-    void initializeGeometry(QVector3D *vertices, int numVertices);
+    void initializeGeometry(QVector3D *vertices);
 
 private:
     QVariant getFromModel(int row, CrystalModel::Columns column) const;
     QMatrix4x4 getCrystalOrientationMatrix() const;
     float getFurthestVertexDistance(QVector3D *vertices, int numVertices) const;
+    void drawAxisLines(QMatrix4x4 viewMatrix, QPainter *painter);
 
     CrystalModel *m_crystals;
     int m_populationIndex;
diff --git a/src/haloray-core/gui/models/simulationStateModel.cpp b/src/haloray-core/gui/models/simulationStateModel.cpp
index 7ca07b4..171f674 100644
--- a/src/haloray-core/gui/models/simulationStateModel.cpp
+++ b/src/haloray-core/gui/models/simulationStateModel.cpp
@@ -36,6 +36,10 @@ SimulationStateModel::SimulationStateModel(SimulationEngine *engine, QObject *pa
     connect(m_simulationEngine, &SimulationEngine::atmosphereChanged, [this]() {
         emit dataChanged(createIndex(0, AtmosphereEnabled), createIndex(0, GroundAlbedo));
     });
+
+    connect(m_simulationEngine, &SimulationEngine::guidesToggled, [this]() {
+        emit dataChanged(createIndex(0, GuidesEnabled), createIndex(0, GuidesEnabled));
+    });
 }
 
 QVariant SimulationStateModel::headerData(int section, Qt::Orientation orientation, int role) const
@@ -74,6 +78,8 @@ QVariant SimulationStateModel::headerData(int section, Qt::Orientation orientati
             return "Atmosphere turbidity";
         case GroundAlbedo:
             return "Ground albedo";
+        case GuidesEnabled:
+            return "Guides enabled";
         }
     }
 
@@ -134,6 +140,8 @@ QVariant SimulationStateModel::data(const QModelIndex &index, int role) const
         return m_simulationEngine->getAtmosphere().turbidity;
     case GroundAlbedo:
         return m_simulationEngine->getAtmosphere().groundAlbedo;
+    case GuidesEnabled:
+        return m_simulationEngine->getGuidesEnabled();
     default:
         break;
     }
@@ -191,6 +199,9 @@ bool SimulationStateModel::setData(const QModelIndex &index, const QVariant &val
     case GroundAlbedo:
         setGroundAlbedo(value.toDouble());
         break;
+    case GuidesEnabled:
+        m_simulationEngine->setGuidesEnabled(value.toBool());
+        break;
     default:
         return false;
     }
diff --git a/src/haloray-core/gui/models/simulationStateModel.h b/src/haloray-core/gui/models/simulationStateModel.h
index e46de67..628be75 100644
--- a/src/haloray-core/gui/models/simulationStateModel.h
+++ b/src/haloray-core/gui/models/simulationStateModel.h
@@ -33,6 +33,7 @@ class SimulationStateModel : public QAbstractTableModel
         AtmosphereEnabled,
         Turbidity,
         GroundAlbedo,
+        GuidesEnabled,
         NUM_COLUMNS
     };
 
diff --git a/src/haloray-core/gui/openGLWidget.cpp b/src/haloray-core/gui/openGLWidget.cpp
index f0f7b33..60c03a6 100644
--- a/src/haloray-core/gui/openGLWidget.cpp
+++ b/src/haloray-core/gui/openGLWidget.cpp
@@ -48,7 +48,7 @@ void OpenGLWidget::paintGL()
     const float adjustedExposure = 500000.0f * m_exposure / (m_engine->getIteration() + 1) / (m_engine->getCamera().fov / 180.0) / m_engine->getRaysPerStep();
     m_textureRenderer->setUniformFloat("adjustedExposure", adjustedExposure);
     m_textureRenderer->setUniformFloat("baseExposure", m_exposure);
-    m_textureRenderer->render(m_engine->getOutputTextureHandle(), m_engine->getBackgroundTextureHandle());
+    m_textureRenderer->render(m_engine->getOutputTextureHandle(), m_engine->getBackgroundTextureHandle(), m_engine->getGuideTextureHandle());
 }
 
 void OpenGLWidget::resizeGL(int w, int h)
diff --git a/src/haloray-core/gui/viewSettingsWidget.cpp b/src/haloray-core/gui/viewSettingsWidget.cpp
index 29d72bb..e3fc356 100644
--- a/src/haloray-core/gui/viewSettingsWidget.cpp
+++ b/src/haloray-core/gui/viewSettingsWidget.cpp
@@ -24,6 +24,7 @@ ViewSettingsWidget::ViewSettingsWidget(SimulationStateModel *viewModel, QWidget
     m_mapper->addMapping(m_pitchSlider, SimulationStateModel::CameraPitch);
     m_mapper->addMapping(m_yawSlider, SimulationStateModel::CameraYaw);
     m_mapper->addMapping(m_hideSubHorizonCheckBox, SimulationStateModel::HideSubHorizon);
+    m_mapper->addMapping(m_showGuidesCheckBox, SimulationStateModel::GuidesEnabled);
     m_mapper->toFirst();
 
     connect(m_cameraProjectionComboBox, QOverload<int>::of(&QComboBox::currentIndexChanged), m_mapper, &QDataWidgetMapper::submit, Qt::QueuedConnection);
@@ -31,6 +32,7 @@ ViewSettingsWidget::ViewSettingsWidget(SimulationStateModel *viewModel, QWidget
     connect(m_pitchSlider, &SliderSpinBox::valueChanged, m_mapper, &QDataWidgetMapper::submit, Qt::QueuedConnection);
     connect(m_yawSlider, &SliderSpinBox::valueChanged, m_mapper, &QDataWidgetMapper::submit, Qt::QueuedConnection);
     connect(m_hideSubHorizonCheckBox, &QCheckBox::stateChanged, m_mapper, &QDataWidgetMapper::submit, Qt::QueuedConnection);
+    connect(m_showGuidesCheckBox, &QCheckBox::stateChanged, m_mapper, &QDataWidgetMapper::submit, Qt::QueuedConnection);
 
     /*
      * It is not possible to map multiple model columns to different properties
@@ -43,7 +45,7 @@ ViewSettingsWidget::ViewSettingsWidget(SimulationStateModel *viewModel, QWidget
     m_maximumFovMapper->toFirst();
 
     connect(m_brightnessSlider, &SliderSpinBox::valueChanged, this, &ViewSettingsWidget::brightnessChanged);
-    connect(m_lockToLightSource, &QCheckBox::stateChanged, this, &ViewSettingsWidget::lockToLightSource);
+    connect(m_lockToLightSourceCheckBox, &QCheckBox::stateChanged, this, &ViewSettingsWidget::lockToLightSource);
 }
 
 void ViewSettingsWidget::setupUi()
@@ -70,7 +72,9 @@ void ViewSettingsWidget::setupUi()
 
     m_hideSubHorizonCheckBox = new QCheckBox();
 
-    m_lockToLightSource = new QCheckBox();
+    m_lockToLightSourceCheckBox = new QCheckBox();
+
+    m_showGuidesCheckBox = new QCheckBox();
 
     auto layout = new QFormLayout(this->contentWidget());
     layout->addRow(tr("Camera projection"), m_cameraProjectionComboBox);
@@ -79,7 +83,8 @@ void ViewSettingsWidget::setupUi()
     layout->addRow(tr("Yaw"), m_yawSlider);
     layout->addRow(tr("Brightness"), m_brightnessSlider);
     layout->addRow(tr("Hide sub-horizon"), m_hideSubHorizonCheckBox);
-    layout->addRow(tr("Lock to light source"), m_lockToLightSource);
+    layout->addRow(tr("Lock to light source"), m_lockToLightSourceCheckBox);
+    layout->addRow(tr("Show guides"), m_showGuidesCheckBox);
 }
 
 void ViewSettingsWidget::setBrightness(double brightness)
diff --git a/src/haloray-core/gui/viewSettingsWidget.h b/src/haloray-core/gui/viewSettingsWidget.h
index 633e836..9ce5fb6 100644
--- a/src/haloray-core/gui/viewSettingsWidget.h
+++ b/src/haloray-core/gui/viewSettingsWidget.h
@@ -32,7 +32,8 @@ class ViewSettingsWidget : public CollapsibleBox
     QComboBox *m_cameraProjectionComboBox;
     QCheckBox *m_hideSubHorizonCheckBox;
     SliderSpinBox *m_brightnessSlider;
-    QCheckBox *m_lockToLightSource;
+    QCheckBox *m_lockToLightSourceCheckBox;
+    QCheckBox *m_showGuidesCheckBox;
 
     SimulationStateModel *m_viewModel;
     QDataWidgetMapper *m_mapper;
diff --git a/src/haloray-core/opengl/texture.cpp b/src/haloray-core/opengl/texture.cpp
index 11ada75..2c9adcf 100644
--- a/src/haloray-core/opengl/texture.cpp
+++ b/src/haloray-core/opengl/texture.cpp
@@ -24,7 +24,7 @@ void Texture::initializeTextureImage()
         glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, m_width, m_height, 0, GL_RGBA, GL_FLOAT, NULL);
         break;
     case Monochrome:
-        glTexImage2D(GL_TEXTURE_2D, 0, GL_R32UI, m_width, m_height, 0, GL_RED, GL_UNSIGNED_INT, NULL);
+        glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, m_width, m_height, 0, GL_RED, GL_FLOAT, NULL);
         break;
     default:
         throw std::runtime_error("Invalid texture type");
diff --git a/src/haloray-core/opengl/textureRenderer.cpp b/src/haloray-core/opengl/textureRenderer.cpp
index 5afc4a9..93967f6 100644
--- a/src/haloray-core/opengl/textureRenderer.cpp
+++ b/src/haloray-core/opengl/textureRenderer.cpp
@@ -71,7 +71,7 @@ void TextureRenderer::initialize()
     m_texDrawProgram = initializeTexDrawShaderProgram();
 }
 
-void TextureRenderer::render(unsigned int haloTextureHandle, int backgroundTextureHandle)
+void TextureRenderer::render(unsigned int haloTextureHandle, unsigned int backgroundTextureHandle, unsigned int guideTextureHandle)
 {
     glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
 
@@ -81,11 +81,18 @@ void TextureRenderer::render(unsigned int haloTextureHandle, int backgroundTextu
     glClearColor(0.0f, 0.0f, 0.0f, 1.0f);
     glClear(GL_COLOR_BUFFER_BIT);
     glBindVertexArray(m_quadVao);
+
     glActiveTexture(GL_TEXTURE0);
     glBindTexture(GL_TEXTURE_2D, haloTextureHandle);
+
     glActiveTexture(GL_TEXTURE1);
     glBindTexture(GL_TEXTURE_2D, backgroundTextureHandle);
+
+    glActiveTexture(GL_TEXTURE2);
+    glBindTexture(GL_TEXTURE_2D, guideTextureHandle);
+
     glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
+
     glActiveTexture(GL_TEXTURE0);
 }
 
diff --git a/src/haloray-core/opengl/textureRenderer.h b/src/haloray-core/opengl/textureRenderer.h
index 6303912..37d5c0a 100644
--- a/src/haloray-core/opengl/textureRenderer.h
+++ b/src/haloray-core/opengl/textureRenderer.h
@@ -13,7 +13,7 @@ class TextureRenderer : protected QOpenGLFunctions_4_4_Core
     void initialize();
     void setUniformFloat(std::string name, float value);
     void render(unsigned int textureHandle);
-    void render(unsigned int haloTextureHandle, int backgroundTextureHandle);
+    void render(unsigned int haloTextureHandle, unsigned int backgroundTextureHandle, unsigned int guideTextureHandle);
     ~TextureRenderer();
 
 private:
diff --git a/src/haloray-core/resources/haloray.qrc b/src/haloray-core/resources/haloray.qrc
index f34cfc5..383b92e 100644
--- a/src/haloray-core/resources/haloray.qrc
+++ b/src/haloray-core/resources/haloray.qrc
@@ -5,6 +5,7 @@
         <file>shaders/sky.glsl</file>
         <file>shaders/renderer.vert</file>
         <file>shaders/renderer.frag</file>
+        <file>shaders/guide.glsl</file>
     </qresource>
     <qresource prefix="/icons">
         <file alias="actions/24/list-add.svg">icons/custom-icons/actions/24/list-add.svg</file>
diff --git a/src/haloray-core/resources/shaders/guide.glsl b/src/haloray-core/resources/shaders/guide.glsl
new file mode 100644
index 0000000..3577072
--- /dev/null
+++ b/src/haloray-core/resources/shaders/guide.glsl
@@ -0,0 +1,122 @@
+#version 440 core
+
+layout(local_size_x = 1, local_size_y = 1) in;
+layout(binding = 2, r32f) uniform image2D outputImage;
+
+uniform struct sunProperties_t
+{
+    float altitude;
+    float diameter;
+} sun;
+
+uniform struct camera_t
+{
+    float pitch;
+    float yaw;
+    float focalLength;
+    int projection;
+    int hideSubHorizon;
+} camera;
+
+#define PROJECTION_STEREOGRAPHIC 0
+#define PROJECTION_RECTILINEAR 1
+#define PROJECTION_EQUIDISTANT 2
+#define PROJECTION_EQUAL_AREA 3
+#define PROJECTION_ORTHOGRAPHIC 4
+
+const float PI = 3.1415926535;
+const float LINEWIDTHDEGREES = 0.25 / sqrt(camera.focalLength);
+const float LINEWIDTH = LINEWIDTHDEGREES * PI / 180.0;
+
+vec2 planarToPolar(vec2 point)
+{
+    float r = length(point);
+    float angle = atan(point.y, point.x);
+    return vec2(r, angle);
+}
+
+mat3 rotateAroundX(float angle)
+{
+    return mat3(
+        1.0, 0.0, 0.0,
+        0.0, cos(angle), sin(angle),
+        0.0, -sin(angle), cos(angle)
+    );
+}
+
+mat3 rotateAroundY(float angle)
+{
+    return mat3(
+        cos(angle), 0.0, -sin(angle),
+        0.0, 1.0, 0.0,
+        sin(angle), 0.0, cos(angle)
+    );
+}
+
+mat3 getCameraOrientationMatrix()
+{
+    return rotateAroundY(-camera.yaw) * rotateAroundX(-camera.pitch);
+}
+
+vec3 getSunVector()
+{
+    return normalize(vec3(0.0, sin(sun.altitude), cos(sun.altitude)));
+}
+
+void main(void)
+{
+    ivec2 resolution = imageSize(outputImage);
+    float aspectRatio = float(resolution.y) / float(resolution.x);
+
+    ivec2 pixelCoordinates = ivec2(gl_GlobalInvocationID.xy);
+
+    vec2 normCoord = pixelCoordinates / vec2(resolution) - 0.5;
+    normCoord.x /= aspectRatio;
+    vec2 polar = planarToPolar(normCoord);
+
+    float projectedAngle;
+
+    // The projection converts 2D coordinates to 3D vectors
+    if (camera.projection == PROJECTION_STEREOGRAPHIC) {
+        projectedAngle = 2.0 * atan(polar.x / 2.0 / camera.focalLength);
+    } else if (camera.projection == PROJECTION_RECTILINEAR) {
+        if (polar.x > 0.5 * PI) return;
+        projectedAngle = atan(polar.x / camera.focalLength);
+    } else if (camera.projection == PROJECTION_EQUIDISTANT) {
+        projectedAngle = polar.x / camera.focalLength;
+    } else if (camera.projection == PROJECTION_EQUAL_AREA) {
+        projectedAngle = 2.0 * asin(polar.x / 2.0 / camera.focalLength);
+    } else if (camera.projection == PROJECTION_ORTHOGRAPHIC) {
+        if (polar.x > 0.5 * PI) return;
+        projectedAngle = asin(polar.x / camera.focalLength);
+    }
+
+    if (projectedAngle > PI) return;
+
+    float x = sin(projectedAngle) * cos(polar.y);
+    float y = sin(projectedAngle) * sin(polar.y);
+    float z = cos(projectedAngle);
+
+    vec3 dir = normalize(getCameraOrientationMatrix() * vec3(x, y, z));
+
+    float horizonDistance = abs(dir.y);
+
+    float zenithDistance = atan(length(dir.xz), dir.y);
+    float nadirDistance = PI - zenithDistance;
+    float minZenithNadirDistance = min(zenithDistance, nadirDistance);
+    const float zenithNadirCrossRadius = 0.025;
+    if (minZenithNadirDistance < zenithNadirCrossRadius) {
+        float minAxisDistance = min(abs(dir.x), abs(dir.z));
+        minZenithNadirDistance = minAxisDistance;
+    }
+
+    float sunDistance = acos(clamp(dot(getSunVector(), dir), -1.0, 1.0));
+    float distance22r = abs(sunDistance - 22.0 / 180.0 * PI);
+    float distance46r = abs(sunDistance - 46.0 / 180.0 * PI);
+    float minRingDistance = min(distance22r, distance46r);
+
+    float minDistance = min(min(minZenithNadirDistance, horizonDistance), minRingDistance);
+    float value = 1.0 - smoothstep(LINEWIDTH * 0.3, LINEWIDTH * 0.5, minDistance);
+
+    imageStore(outputImage, pixelCoordinates, vec4(value, 0.0, 0.0, 0.0));
+}
diff --git a/src/haloray-core/resources/shaders/raytrace.glsl b/src/haloray-core/resources/shaders/raytrace.glsl
index 1d63cb8..9a5fa17 100644
--- a/src/haloray-core/resources/shaders/raytrace.glsl
+++ b/src/haloray-core/resources/shaders/raytrace.glsl
@@ -79,7 +79,7 @@ ivec3 triangles[] = ivec3[](
     ivec3(0, 3, 4),
     ivec3(0, 4, 5),
 
-    // Face 1 (pyramid edges)
+    // Upper pyramid faces
     ivec3(0, 6, 1),
     ivec3(6, 7, 1),
     ivec3(1, 7, 2),
@@ -99,7 +99,7 @@ ivec3 triangles[] = ivec3[](
     ivec3(18, 22, 21),
     ivec3(18, 23, 22),
 
-    // Face 2 (pyramid edges)
+    // Lower pyramid faces
     ivec3(12, 18, 13),
     ivec3(18, 19, 13),
     ivec3(13, 19, 14),
@@ -140,6 +140,10 @@ ivec3 triangles[] = ivec3[](
 
 vec3 triangleNormalCache[triangles.length()];
 
+// ***********************************
+// Random number generator functions *
+// ***********************************
+
 uint wang_hash(uint a)
 {
     a -= (a << 6);
@@ -172,6 +176,10 @@ vec2 randn(void)
     return vec2(u1 * cos(u2), u1 * sin(u2));
 }
 
+// **************************
+// Color matching functions *
+// **************************
+
 float xFit_1931(float wave)
 {
     float t1 = (wave - 442.0) * ((wave < 442.0) ? 0.0624 : 0.0374);
@@ -194,6 +202,10 @@ float zFit_1931(float wave)
     return 1.217 * exp(-0.5 * t1 * t1) + 0.681 * exp(-0.5 * t2 * t2);
 }
 
+// ***********************
+// Ray tracing functions *
+// ***********************
+
 float getIceIOR(float wavelength)
 {
     // Eq. from Simulating rainbows and halos in color by Stanley Gedzelman
@@ -211,7 +223,7 @@ uint selectFirstTriangle(vec3 rayDirection)
         vec3 v0 = vertices[triangle.x];
         vec3 v1 = vertices[triangle.y];
         vec3 v2 = vertices[triangle.z];
-        vec3 triangleCrossProduct = cross(v2 - v0, v1 - v0);
+        vec3 triangleCrossProduct = cross(v1 - v0, v2 - v0);
         float triangleArea = 0.5 * length(triangleCrossProduct);
         vec3 triangleNormal = normalize(triangleCrossProduct);
         triangleNormalCache[i] = triangleNormal;
@@ -252,7 +264,7 @@ vec3 sampleTriangle(uint triangleIndex)
 
 vec3 getNormal(uint triangleIndex)
 {
-    return -triangleNormalCache[triangleIndex];
+    return triangleNormalCache[triangleIndex];
 }
 
 float getReflectionCoefficient(vec3 normal, vec3 rayDir, float n0, float n1)
@@ -269,14 +281,22 @@ float getReflectionCoefficient(vec3 normal, vec3 rayDir, float n0, float n1)
     return 0.5 * (rs + rp);
 }
 
+// Find triangle which ray intersects with using
+// the Trumbore-Möller algorithm.
 intersection findIntersection(vec3 rayOrigin, vec3 rayDirection)
 {
     for (int triangleIndex = 0; triangleIndex < triangles.length(); ++triangleIndex)
     {
+        // Note the ordering of the vertices!
+        // v1 and v2 are switched, because with
+        // the conventional ordering the triangle
+        // normal ends up pointing outwards, and it
+        // must point inwards when tracing rays
+        // inside the ice crystal.
         ivec3 triangle = triangles[triangleIndex];
         vec3 v0 = vertices[triangle.x];
-        vec3 v1 = vertices[triangle.y];
-        vec3 v2 = vertices[triangle.z];
+        vec3 v1 = vertices[triangle.z];
+        vec3 v2 = vertices[triangle.y];
 
         vec3 v0v1 = v1 - v0;
         vec3 v0v2 = v2 - v0;
@@ -309,7 +329,7 @@ vec3 traceRay(vec3 rayOrigin, vec3 rayDirection, float indexOfRefraction)
     {
         intersection hitResult = findIntersection(ro, rd);
         if (hitResult.didHit == false) break;
-        vec3 normal = getNormal(hitResult.triangleIndex);
+        vec3 normal = -getNormal(hitResult.triangleIndex);
         float reflectionCoefficient = getReflectionCoefficient(normal, rd, indexOfRefraction, 1.0);
         if (rand() < reflectionCoefficient)
         {
@@ -324,6 +344,31 @@ vec3 traceRay(vec3 rayOrigin, vec3 rayDirection, float indexOfRefraction)
     return vec3(0.0);
 }
 
+vec3 castRayThroughCrystal(vec3 rayDirection, float wavelength)
+{
+    uint triangleIndex = selectFirstTriangle(rayDirection);
+    vec3 startingPoint = sampleTriangle(triangleIndex);
+    vec3 startingPointNormal = getNormal(triangleIndex);
+    float indexOfRefraction = getIceIOR(wavelength);
+    float reflectionCoeff = getReflectionCoefficient(startingPointNormal, rayDirection, 1.0, indexOfRefraction);
+    vec3 resultRay = vec3(0.0);
+    if (rand() < reflectionCoeff)
+    {
+        // Ray reflects off crystal
+        resultRay = reflect(rayDirection, startingPointNormal);
+    } else {
+        // Ray enters crystal
+        vec3 refractedRayDirection = refract(rayDirection, startingPointNormal, 1.0 / indexOfRefraction);
+        resultRay = traceRay(startingPoint, refractedRayDirection, indexOfRefraction);
+    }
+
+    return resultRay;
+}
+
+// **************************************
+// Sun direction and spectrum functions *
+// **************************************
+
 vec3 getSunDirection(float altitude)
 {
     // X and Z are horizontal, sun moves on the Y-Z plane
@@ -349,6 +394,22 @@ vec3 sampleSun(float altitude)
     return normalize(sampleDirection);
 }
 
+float daylightEstimate(float wavelength)
+{
+    return 1.0 - 0.0013333 * wavelength;
+}
+
+float sampleSunSpectrum(float wavelength)
+{
+    int index = clamp(int(floor((wavelength - 400.0) / 10.0)), 0, 29);
+    float wavelengthFract = (wavelength - (400.0 + index * 10.0)) / 10.0;
+    return mix(sun.spectrum[index], sun.spectrum[index + 1], wavelengthFract);
+}
+
+// ********************
+// Rotation functions *
+// ********************
+
 mat3 rotateAroundX(float angle)
 {
     return mat3(
@@ -376,9 +437,18 @@ mat3 rotateAroundZ(float angle)
     );
 }
 
+vec2 rotate2D(float angle, vec2 point)
+{
+    return mat2(cos(angle), sin(angle), -sin(angle), cos(angle)) * point;
+}
+
+// ******************************************
+// Camera and crystal orientation functions *
+// ******************************************
+
 mat3 getCameraOrientationMatrix()
 {
-    return rotateAroundX(camera.pitch) * rotateAroundY(camera.yaw);
+    return rotateAroundY(-camera.yaw) * rotateAroundX(-camera.pitch);
 }
 
 mat3 getUniformRandomRotationMatrix(void)
@@ -392,13 +462,6 @@ mat3 getUniformRandomRotationMatrix(void)
     return (2.0 * outerProduct(reflectionVector, reflectionVector) - mat3(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) * zRotationMatrix;
 }
 
-vec2 cartesianToPolar(vec3 direction)
-{
-    float r = atan(length(direction.xy), direction.z);
-    float angle = atan(direction.y, direction.x);
-    return vec2(r, angle);
-}
-
 mat3 getRotationMatrix(void)
 {
     if (crystalProperties.tiltDistribution == DISTRIBUTION_UNIFORM && crystalProperties.rotationDistribution == DISTRIBUTION_UNIFORM)
@@ -409,18 +472,18 @@ mat3 getRotationMatrix(void)
     // Tilt of the crystal C-axis
     mat3 tiltMat;
 
-    // Rotation around crystal C-axis
-    mat3 rotationMat;
-
     if (crystalProperties.tiltDistribution == DISTRIBUTION_UNIFORM) {
         tiltMat = rotateAroundZ(rand() * 2.0 * PI);
     } else {
         float angleAverage = crystalProperties.tiltAverage;
         float angleStd = crystalProperties.tiltStd;
         float tiltAngle = angleAverage + angleStd * randn().x;
-        tiltMat = rotateAroundZ(tiltAngle);
+        tiltMat = rotateAroundZ(-tiltAngle);
     }
 
+    // Rotation around crystal C-axis
+    mat3 rotationMat;
+
     if (crystalProperties.rotationDistribution == DISTRIBUTION_UNIFORM)
     {
         rotationMat = rotateAroundY(rand() * 2.0 * PI);
@@ -434,224 +497,280 @@ mat3 getRotationMatrix(void)
     return rotateAroundY(rand() * 2.0 * PI) * tiltMat * rotationMat;
 }
 
-float daylightEstimate(float wavelength)
-{
-    return 1.0 - 0.0013333 * wavelength;
-}
+// ***************************************
+// Crystal geometry generation functions *
+// ***************************************
 
-float sampleSunSpectrum(float wavelength)
+int getNext(int i)
 {
-    int index = clamp(int(floor((wavelength - 400.0) / 10.0)), 0, 29);
-    float wavelengthFract = (wavelength - (400.0 + index * 10.0)) / 10.0;
-    return mix(sun.spectrum[index], sun.spectrum[index + 1], wavelengthFract);
+    return int(mod(i + 1, 6));
 }
 
-void storePixel(ivec2 pixelCoordinates, vec3 value)
+int getPrev(int i)
 {
-    memoryBarrierImage();
-    vec3 currentValue = imageLoad(outputImage, pixelCoordinates).xyz;
-    imageStore(outputImage, pixelCoordinates, vec4(currentValue + value, 1.0));
+    return int(mod(mod(i - 1, 6) + 6, 6));
 }
 
-vec3 castRayThroughCrystal(vec3 rayDirection, float wavelength)
+vec3[6] generateApexNormals(float apexAngle)
 {
-    uint triangleIndex = selectFirstTriangle(rayDirection);
-    vec3 startingPoint = sampleTriangle(triangleIndex);
-    vec3 startingPointNormal = -getNormal(triangleIndex);
-    float indexOfRefraction = getIceIOR(wavelength);
-    float reflectionCoeff = getReflectionCoefficient(startingPointNormal, rayDirection, 1.0, indexOfRefraction);
-    vec3 resultRay = vec3(0.0);
-    if (rand() < reflectionCoeff)
+    vec3 apexNormals[6];
+    for (int i = 0; i < 6; ++i)
     {
-        // Ray reflects off crystal
-        resultRay = reflect(rayDirection, startingPointNormal);
-    } else {
-        // Ray enters crystal
-        vec3 refractedRayDirection = refract(rayDirection, startingPointNormal, 1.0 / indexOfRefraction);
-        resultRay = traceRay(startingPoint, refractedRayDirection, indexOfRefraction);
+        float rotAngle = i * PI / 3.0;
+        float halfApex = apexAngle / 2.0;
+        apexNormals[i] = rotateAroundY(rotAngle) * rotateAroundX(-halfApex) * vec3(0.0, 0.0, 1.0);
     }
-
-    return resultRay;
+    return apexNormals;
 }
 
-/* Lines are represented in Hesse normal form, where X component
-   of the vector is the closest distance from origin to the line,
-   and Y component of the vector is the angle of the line's normal
-   in radians. */
-vec2 lineIntersect(vec2 line1, vec2 line2)
+float getMaximumApexHeight(vec3 normals[6],
+                           vec3 vertices[24],
+                           float prismFaceDistances[6],
+                           float apexAngle,
+                           int vertexOffset)
 {
-    float p1 = line1.x;
-    float theta1 = line1.y;
-
-    float p2 = line2.x;
-    float theta2 = line2.y;
+    float maxApexHeight = 1e38;
+    // Find non-convex pyramid face heights
+    for (int face = 0; face < 6; ++face)
+    {
+        int prevFace = getPrev(face);
+        int nextFace = getNext(face);
+
+        vec3 nPrev = normals[prevFace];
+        vec3 nCurr = normals[face];
+        vec3 nNext = normals[nextFace];
+
+        vec3 prevVert = vertices[prevFace + vertexOffset];
+        vec3 nextVert = vertices[face + vertexOffset];
+
+        vec3 numerator = dot(prevVert, nPrev) * cross(nCurr, nNext) +
+                         dot(prevVert, nCurr) * cross(nNext, nPrev) +
+                         dot(nextVert, nNext) * cross(nPrev, nCurr);
+        mat3 detMatrix = mat3(nPrev, nCurr, nNext);
+        float denominator = determinant(detMatrix);
+        float faceHeight = abs((numerator / denominator).y);
+        maxApexHeight = min(faceHeight, maxApexHeight);
+    }
 
-    float deltaSine = sin(theta2 - theta1);
-    float x = (p1 * sin(theta2) - p2 * sin(theta1)) / deltaSine;
-    float y = (p2 * cos(theta1) - p1 * cos(theta2)) / deltaSine;
+    // Find maximum apex heights based on prism face distances and apex angle
+    for (int i = 0; i < 3; ++i)
+    {
+        float dist = prismFaceDistances[i];
+        float distOpposite = prismFaceDistances[i + 3];
+        float h = (dist + distOpposite) / (2.0 * tan(apexAngle / 2.0));
+        maxApexHeight = min(h, maxApexHeight);
+    }
 
-    return vec2(x, y);
+    return maxApexHeight;
 }
 
 void initializeCrystal()
 {
-    float deltaAngle = radians(60.0);
     vec2 hexagonCorners[6];
-    /* The sqrt(3)/2 multiplier makes the default crystal such
-       that the distance of a vertex from the C axis is 1.0.
-       sqrt(3)/2 equals cos(30deg), which is faster to calculate
-       on GPU */
-    float sizeScaler = cos(radians(30.0));
-    for (int face = 0; face < 6; ++face)
+
+    // Calculate initial hexagon shape
+    for (int i = 0; i < 6; ++i)
     {
-        int previousFace = face == 0 ? 5 : face - 1;
-        int nextFace = face == 5 ? 0 : face + 1;
+        float d1 = crystalProperties.prismFaceDistances[i];
+        float d2 = crystalProperties.prismFaceDistances[getNext(i)];
 
-        float previousAngle = (face + 1) * deltaAngle;
-        float currentAngle = previousAngle + deltaAngle;
-        float nextAngle = previousAngle + 2.0 * deltaAngle;
+        float angle = -i * PI / 3.0;
+        float x_stat = 2.0 * d2 / sqrt(3.0) - d1 / sqrt(3.0);
+        float y_stat = d1;
+        hexagonCorners[i] = rotate2D(angle, vec2(x_stat, y_stat));
+    }
 
-        float previousDistance = sizeScaler * crystalProperties.prismFaceDistances[previousFace];
-        float currentDistance = sizeScaler * crystalProperties.prismFaceDistances[face];
-        float nextDistance = sizeScaler * crystalProperties.prismFaceDistances[nextFace];
+    // Fix degenerate edges on hexagon
+    for (int face = 0; face < 6; ++face)
+    {
+        int prevFace = getPrev(face);
+        int nextFace = getNext(face);
+        float angle = -face * PI / 3.0;
+
+        float d1 = crystalProperties.prismFaceDistances[face];
+        float d2 = crystalProperties.prismFaceDistances[nextFace];
+        float d3 = crystalProperties.prismFaceDistances[prevFace];
+        if (d1 > d2 + d3)
+        {
+            float x_stat = d2 / sqrt(3.0) - d3 / sqrt(3.0);
+            float y_stat = d2 + d3;
+            vec2 rotatedPoint = rotate2D(angle, vec2(x_stat, y_stat));
+            hexagonCorners[face] = rotatedPoint;
+            hexagonCorners[prevFace] = rotatedPoint;
+        }
+    }
 
-        vec2 previousLine = vec2(previousDistance, previousAngle);
-        vec2 currentLine = vec2(currentDistance, currentAngle);
-        vec2 nextLine = vec2(nextDistance, nextAngle);
+    /* Scaling value makes sure eventual A axis length is 2.0, so C/A ratio
+     * can be easily corrected. */
+    float hexagonScaler = length(hexagonCorners[1] - hexagonCorners[4]);
+    for (int i = 0; i < 6; ++i)
+    {
+        hexagonCorners[i] *= 2.0 / hexagonScaler;
+    }
 
-        vec2 previousCurrentIntersection = lineIntersect(previousLine, currentLine);
-        vec2 currentNextIntersection = lineIntersect(currentLine, nextLine);
-        vec2 previousNextIntersection = lineIntersect(previousLine, nextLine);
+    for (int face = 0; face < 6; ++face)
+    {
+        // Corresponding vertices of each crystal layer have the same X and Z coordinates
+        // First six vertices are the top apex cap
+        // Next six vertices are the top of the base hexagonal crystal
+        // Next six vertices are the bottom of the base hexagonal crystal
+        // Last six vertices are the bottom apex cap
+        vec3 v = vec3(hexagonCorners[face].x, 0.0, hexagonCorners[face].y);
+        vertices[face] = v;
+        vertices[face + 6] = v;
+        vertices[face + 12] = v;
+        vertices[face + 18] = v;
+    }
 
-        float previousCurrentIntersectionDistance = length(previousCurrentIntersection);
-        float currentNextIntersectionDistance = length(currentNextIntersection);
-        float previousNextIntersectionDistance = length(previousNextIntersection);
+    vec2 random = randn();
+    float upperApexHeight = clamp(crystalProperties.upperApexHeightAverage + crystalProperties.upperApexHeightStd * random.x, 0.0, 1.0);
+    float lowerApexHeight = clamp(crystalProperties.lowerApexHeightAverage + crystalProperties.lowerApexHeightStd * random.y, 0.0, 1.0);
 
-        vec2 v1 = previousCurrentIntersectionDistance < previousNextIntersectionDistance ? previousCurrentIntersection : previousNextIntersection;
-        vec2 v2 = currentNextIntersectionDistance < previousNextIntersectionDistance ? currentNextIntersection : previousNextIntersection;
+    if (upperApexHeight > 0.0 && crystalProperties.upperApexAngle < PI && crystalProperties.upperApexAngle > 0.0)
+    {
+        // Generate normals for upper pyramid cap
+        vec3 upperApexNormals[6] = generateApexNormals(crystalProperties.upperApexAngle);
 
-        if (face > 0 && previousNextIntersectionDistance > length(hexagonCorners[face]))
+        // Set upper pyramid cap vertex positions
+        float maxUpperApexHeight = getMaximumApexHeight(upperApexNormals, vertices, crystalProperties.prismFaceDistances, crystalProperties.upperApexAngle, 0);
+        for (int i = 0; i < 6; ++i)
         {
-            v1 = hexagonCorners[face];
+            int next = getNext(i);
+            vec3 pyramidEdge = cross(upperApexNormals[i], upperApexNormals[next]);
+            vertices[i] += upperApexHeight * maxUpperApexHeight * pyramidEdge / pyramidEdge.y;
         }
+    }
 
-        if (face == 5 && previousNextIntersectionDistance > length(hexagonCorners[nextFace]))
+    if (lowerApexHeight > 0.0 && crystalProperties.lowerApexAngle < PI && crystalProperties.lowerApexAngle > 0.0)
+    {
+        // Generate normals for lower pyramid cap
+        vec3 lowerApexNormals[6] = generateApexNormals(crystalProperties.lowerApexAngle);
+        for (int i = 0; i < 6; ++i)
         {
-            v2 = hexagonCorners[nextFace];
+            lowerApexNormals[i].y *= -1.0;
         }
 
-        hexagonCorners[face] = v1;
-        hexagonCorners[nextFace] = v2;
+        // Set lower pyramid cap vertex positions
+        float maxLowerApexHeight = getMaximumApexHeight(lowerApexNormals, vertices, crystalProperties.prismFaceDistances, crystalProperties.lowerApexAngle, 18);
+        for (int i = 0; i < 6; ++i)
+        {
+            int next = getNext(i);
+            vec3 pyramidEdge = cross(lowerApexNormals[i], lowerApexNormals[next]);
+            vertices[i + 18] -= lowerApexHeight * maxLowerApexHeight * pyramidEdge / pyramidEdge.y;
+        }
     }
 
-    for (int face = 0; face < 6; ++face)
+    // Scale crystal vertically to have correct C/A ratio
+    float caRatio = max(0.0, crystalProperties.caRatioAverage + randn().x * crystalProperties.caRatioStd);
+    for (int i = 0; i < 12; ++i)
     {
-        vertices[face].xz = hexagonCorners[face];
-        vertices[face].y = 1.0;
-
-        vertices[face + 6].xz = hexagonCorners[face];
-        vertices[face + 6].y = 1.0;
-
-        vertices[face + 12].xz = hexagonCorners[face];
-        vertices[face + 12].y = -1.0;
-
-        vertices[face + 18].xz = hexagonCorners[face];
-        vertices[face + 18].y = -1.0;
+        vertices[i].y += caRatio;
+        vertices[i + 12].y -= caRatio;
     }
 
-    // Stretch the crystal to correct C/A ratio
-    float caMultiplier = max(0.0, crystalProperties.caRatioAverage + randn().x * crystalProperties.caRatioStd);
-    for (int i = 0; i < vertices.length(); ++i)
+    // Rotate crystal around C-axis so that face numbering follows conventions
+    // Prism face 0 (Face 3 in the UI) should be up in a column Parry position
+    mat3 conventionMatrix = rotateAroundY(-PI / 2.0);
+    for (int i = 0; i < 24; ++i)
     {
-        vertices[i].y *= max(0.0, caMultiplier);
+        vertices[i] = conventionMatrix * vertices[i];
     }
+}
 
-    // Scale pyramid caps
-    float upperApexMaxHeight = sizeScaler / tan(crystalProperties.upperApexAngle / 2.0);
-    float lowerApexMaxHeight = sizeScaler / tan(crystalProperties.lowerApexAngle / 2.0);
-
-    vec2 random = randn();
-    float upperApexHeight = clamp(crystalProperties.upperApexHeightAverage + crystalProperties.upperApexHeightStd * random.x, 0.0, 1.0);
-    float lowerApexHeight = clamp(crystalProperties.lowerApexHeightAverage + crystalProperties.lowerApexHeightStd * random.y, 0.0, 1.0);
+// ***********
+// Utilities *
+// ***********
 
-    for (int i = 0; i < 6; ++i)
-    {
-        vertices[i].xz *= 1.0 - upperApexHeight;
-        vertices[i].y += upperApexHeight * upperApexMaxHeight;
+// Convert 3D unit vector to polar coordinates in the XY plane
+vec2 cartesianToPolar(vec3 direction)
+{
+    float r = atan(length(direction.xy), direction.z);
+    float angle = atan(direction.y, direction.x);
+    return vec2(r, angle);
+}
 
-        vertices[vertices.length() - i - 1].xz *= 1.0 - lowerApexHeight;
-        vertices[vertices.length() - i - 1].y -= lowerApexHeight * lowerApexMaxHeight;
-    }
+void storePixel(ivec2 pixelCoordinates, vec3 value)
+{
+    memoryBarrierImage();
+    vec3 currentValue = imageLoad(outputImage, pixelCoordinates).xyz;
+    imageStore(outputImage, pixelCoordinates, vec4(currentValue + value, 1.0));
 }
 
+// ************
+// Main logic *
+// ************
+
 void main(void)
 {
     initializeCrystal();
 
-    vec3 rayDirection = -sampleSun(sun.altitude);
+    // Pick ray wavelength in nanometers
     float wavelength = 400.0 + rand() * 300.0;
 
-    // Rotation matrix to orient ray/crystal
-    mat3 rotationMatrix = getRotationMatrix();
+    // Generate incoming ray from sun (points from sun to origin)
+    vec3 incidentSunRay = -sampleSun(sun.altitude);
 
-    /* The inverse rotation matrix must be applied because we are
-    rotating the incoming ray and not the crystal itself. */
-    vec3 rotatedRayDirection = normalize(rayDirection * rotationMatrix);
-
-    vec3 resultRay = castRayThroughCrystal(rotatedRayDirection, wavelength);
+    // The inverse rotation matrix must be applied because we are
+    // rotating the incoming ray and not the crystal itself. The ray is
+    // now transformed into crystal-oriented space.
+    mat3 rotationMatrix = getRotationMatrix();
+    vec3 incidentRayCrystalSpace = normalize(transpose(rotationMatrix) * incidentSunRay);
 
-    if (length(resultRay) < 0.0001) return;
+    vec3 exitantRayCrystalSpace = castRayThroughCrystal(incidentRayCrystalSpace, wavelength);
+    if (length(exitantRayCrystalSpace) < 0.0001) return;
 
-    resultRay = rotationMatrix * resultRay;
+    // Return ray back to world coordinate space
+    vec3 exitantRay = rotationMatrix * exitantRayCrystalSpace;
 
     if (multipleScatter != 0.0 && multipleScatter > rand())
     {
-        // Rotation matrix to orient ray/crystal
+        // The inverse rotation matrix must be applied because we are
+        // rotating the incoming ray and not the crystal itself. The ray is
+        // now transformed into crystal-oriented space.
         rotationMatrix = getRotationMatrix();
+        vec3 secondIncidentRayCrystalSpace = normalize(transpose(rotationMatrix) * exitantRay);
 
-        /* The inverse rotation matrix must be applied because we are
-        rotating the incoming ray and not the crystal itself. */
-        rotatedRayDirection = normalize(resultRay * rotationMatrix);
+        vec3 secondExitantRayCrystalSpace = castRayThroughCrystal(secondIncidentRayCrystalSpace, wavelength);
+        if (length(secondExitantRayCrystalSpace) < 0.0001) return;
 
-        resultRay = castRayThroughCrystal(rotatedRayDirection, wavelength);
-
-        if (length(resultRay) < 0.0001) return;
-
-        resultRay = rotationMatrix * resultRay;
+        // Return ray back to world coordinate space
+        exitantRay = rotationMatrix * secondExitantRayCrystalSpace;
     }
 
     // Hide subhorizon rays
-    if (camera.hideSubHorizon == 1 && resultRay.y > 0.0) return;
+    if (camera.hideSubHorizon == 1 && exitantRay.y > 0.0) return;
 
     ivec2 resolution = imageSize(outputImage);
     float aspectRatio = float(resolution.y) / float(resolution.x);
 
-    resultRay = normalize(-getCameraOrientationMatrix() * resultRay);
-    vec2 polar = cartesianToPolar(resultRay);
+    // Camera is looking down the positive Z axis.
+    // Ray is now transformed into camera space.
+    // Camera matrix must be inverted (transposed) because ray
+    // is being transformed, not the camera.
+    vec3 exitantRayCameraSpace = normalize(transpose(getCameraOrientationMatrix()) * exitantRay);
+    vec3 lightDirectionCameraSpace = -exitantRayCameraSpace;
+    vec2 polar = cartesianToPolar(lightDirectionCameraSpace);
+    float polarRadius = polar.x;
+    float polarAngle = polar.y;
 
     float projectionFunction;
-
-    // The projection converts 3D vectors to 2D points
     if (camera.projection == PROJECTION_STEREOGRAPHIC) {
-        projectionFunction = 2.0 * tan(polar.x / 2.0);
+        projectionFunction = 2.0 * tan(polarRadius / 2.0);
     } else if (camera.projection == PROJECTION_RECTILINEAR) {
-        if (polar.x > 0.5 * PI) return;
-        projectionFunction = tan(polar.x);
+        if (polarRadius > 0.5 * PI) return;
+        projectionFunction = tan(polarRadius);
     } else if (camera.projection == PROJECTION_EQUIDISTANT) {
-        projectionFunction = polar.x;
+        projectionFunction = polarRadius;
     } else if (camera.projection == PROJECTION_EQUAL_AREA) {
-        projectionFunction = 2.0 * sin(polar.x / 2.0);
+        projectionFunction = 2.0 * sin(polarRadius / 2.0);
     } else if (camera.projection == PROJECTION_ORTHOGRAPHIC) {
-        if (polar.x > 0.5 * PI) return;
-        projectionFunction = sin(polar.x);
+        if (polarRadius > 0.5 * PI) return;
+        projectionFunction = sin(polarRadius);
     }
 
-    vec2 projected = camera.focalLength * projectionFunction * vec2(aspectRatio * cos(polar.y), sin(polar.y));
+    vec2 projected = camera.focalLength * projectionFunction * vec2(aspectRatio * cos(polarAngle), sin(polarAngle));
     vec2 normalizedCoordinates = 0.5 + projected;
 
-    if (any(lessThanEqual(normalizedCoordinates, vec2(0.0))) || any(greaterThanEqual(normalizedCoordinates, vec2(1.0))))
-        return;
-
     float sunRadiance;
     if (atmosphereEnabled == 1)
     {
diff --git a/src/haloray-core/resources/shaders/renderer.frag b/src/haloray-core/resources/shaders/renderer.frag
index 7e5db17..82090d7 100644
--- a/src/haloray-core/resources/shaders/renderer.frag
+++ b/src/haloray-core/resources/shaders/renderer.frag
@@ -5,6 +5,7 @@ uniform float baseExposure;
 uniform float adjustedExposure;
 layout (binding = 0) uniform sampler2D haloTexture;
 layout (binding = 1) uniform sampler2D backgroundTexture;
+layout (binding = 2) uniform sampler2D guideTexture;
 
 #define FXAA_REDUCE_MIN   (1.0/ 128.0)
 #define FXAA_REDUCE_MUL   (1.0 / 8.0)
@@ -69,7 +70,11 @@ void main(void) {
     vec4 antialiasedBackground = fxaa(backgroundTexture, gl_FragCoord.xy);
     vec3 backgroundLinearSrgb = max(vec3(0.0), baseExposure * antialiasedBackground.rgb);
     vec3 haloLinearSrgb = adjustedExposure * texelFetch(haloTexture, ivec2(gl_FragCoord.xy), 0).xyz;
+
     vec3 linearImage = 0.005 * backgroundLinearSrgb + 0.1 * haloLinearSrgb;
     vec3 gammaCorrected = 1.055 * pow(linearImage, vec3(0.417)) - 0.055;
-    color = vec4(clamp(gammaCorrected, 0.0, 1.0), 1.0);
+    vec3 guide = vec3(texelFetch(guideTexture, ivec2(gl_FragCoord.xy), 0).r);
+
+    vec3 finalColor = clamp(gammaCorrected + 0.5 * guide, 0.0, 1.0);
+    color = vec4(finalColor, 1.0);
 }
diff --git a/src/haloray-core/resources/shaders/sky.glsl b/src/haloray-core/resources/shaders/sky.glsl
index fa0c7d4..e99ab0a 100644
--- a/src/haloray-core/resources/shaders/sky.glsl
+++ b/src/haloray-core/resources/shaders/sky.glsl
@@ -1,7 +1,7 @@
 #version 440 core
 
 layout(local_size_x = 1, local_size_y = 1) in;
-layout(binding = 2, rgba32f) uniform coherent image2D outputImage;
+layout(binding = 1, rgba32f) uniform image2D outputImage;
 
 uniform struct sunProperties_t
 {
@@ -43,10 +43,7 @@ const float mixingMinElevation = radians(0.0);
 
 vec3 getSunVector()
 {
-    /* NOTE: The sun vector is now in the opposite Z direction
-      than in the crystal raytracing shader. This should probably
-      made the same in all shaders. */
-    return normalize(vec3(0.0, sin(sun.altitude), -cos(sun.altitude)));
+    return normalize(vec3(0.0, sin(sun.altitude), cos(sun.altitude)));
 }
 
 /*
@@ -230,7 +227,7 @@ mat3 rotateAroundY(float angle)
 
 mat3 getCameraOrientationMatrix()
 {
-    return rotateAroundY(camera.yaw) * rotateAroundX(camera.pitch);
+    return rotateAroundY(-camera.yaw) * rotateAroundX(-camera.pitch);
 }
 
 vec3 renderSun(vec3 direction)
@@ -290,7 +287,7 @@ void main(void)
 
     float x = sin(projectedAngle) * cos(polar.y);
     float y = sin(projectedAngle) * sin(polar.y);
-    float z = -cos(projectedAngle);
+    float z = cos(projectedAngle);
 
     vec3 dir = normalize(getCameraOrientationMatrix() * vec3(x, y, z));
     if (dir.y < 0.0) return;
diff --git a/src/haloray-core/simulation/simulationEngine.cpp b/src/haloray-core/simulation/simulationEngine.cpp
index b34c00c..05869f6 100644
--- a/src/haloray-core/simulation/simulationEngine.cpp
+++ b/src/haloray-core/simulation/simulationEngine.cpp
@@ -33,7 +33,8 @@ SimulationEngine::SimulationEngine(
       m_cameraLockedToLightSource(false),
       m_multipleScatteringProbability(0.0),
       m_crystalRepository(crystalRepository),
-      m_atmosphere(Atmosphere::createDefaultAtmosphere())
+      m_atmosphere(Atmosphere::createDefaultAtmosphere()),
+      m_guidesEnabled(false)
 {
     initialize();
 }
@@ -94,6 +95,20 @@ void SimulationEngine::setAtmosphere(Atmosphere atmosphere)
     emit atmosphereChanged(m_atmosphere);
 }
 
+bool SimulationEngine::getGuidesEnabled() const
+{
+    return m_guidesEnabled;
+}
+
+void SimulationEngine::setGuidesEnabled(bool newState)
+{
+    if (m_guidesEnabled == newState) return;
+
+    clear();
+    m_guidesEnabled = newState;
+    emit guidesToggled(m_guidesEnabled);
+}
+
 unsigned int SimulationEngine::getOutputTextureHandle() const
 {
     return m_simulationTexture->getHandle();
@@ -104,6 +119,11 @@ unsigned int SimulationEngine::getBackgroundTextureHandle() const
     return m_backgroundTexture->getHandle();
 }
 
+unsigned int SimulationEngine::getGuideTextureHandle() const
+{
+    return m_guideTexture->getHandle();
+}
+
 unsigned int SimulationEngine::getIteration() const
 {
     return m_iteration;
@@ -127,6 +147,24 @@ void SimulationEngine::step()
 {
     ++m_iteration;
 
+    if (m_guidesEnabled && m_iteration == 1)
+    {
+        glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
+        glBindImageTexture(m_guideTexture->getTextureUnit(), m_guideTexture->getHandle(), 0, GL_FALSE, 0, GL_READ_WRITE, GL_R32F);
+        m_guideShader->bind();
+        glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
+
+        m_guideShader->setUniformValue("sun.altitude", degToRad(m_light.altitude));
+        m_guideShader->setUniformValue("sun.diameter", degToRad(m_light.diameter));
+        m_guideShader->setUniformValue("camera.pitch", degToRad(m_camera.pitch));
+        m_guideShader->setUniformValue("camera.yaw", degToRad(m_camera.yaw));
+        m_guideShader->setUniformValue("camera.focalLength", m_camera.getFocalLength());
+        m_guideShader->setUniformValue("camera.projection", m_camera.projection);
+        m_guideShader->setUniformValue("camera.hideSubHorizon", m_camera.hideSubHorizon ? 1 : 0);
+
+        glDispatchCompute(m_outputWidth, m_outputHeight, 1);
+    }
+
     if (m_atmosphere.enabled && m_iteration == 1)
     {
         auto skyState = SkyModel::Create(degToRad(m_light.altitude), m_atmosphere.turbidity, m_atmosphere.groundAlbedo, degToRad(m_light.diameter / 2.0));
@@ -231,6 +269,10 @@ void SimulationEngine::clear()
 
     glClearTexImage(m_backgroundTexture->getHandle(), 0, GL_RGBA, GL_FLOAT, NULL);
     glBindImageTexture(m_backgroundTexture->getTextureUnit(), m_backgroundTexture->getHandle(), 0, GL_FALSE, 0, GL_READ_WRITE, GL_RGBA32F);
+
+    glClearTexImage(m_guideTexture->getHandle(), 0, GL_RGBA, GL_FLOAT, NULL);
+    glBindImageTexture(m_guideTexture->getTextureUnit(), m_guideTexture->getHandle(), 0, GL_FALSE, 0, GL_READ_WRITE, GL_R32F);
+
     m_iteration = 0;
 }
 
@@ -262,7 +304,7 @@ void SimulationEngine::initialize()
 void SimulationEngine::initializeShaders()
 {
     qInfo("Initializing raytracing shader");
-    m_simulationShader = std::make_unique<QOpenGLShaderProgram>();
+    m_simulationShader = new QOpenGLShaderProgram(this);
     bool raytraceShaderReadSucceeded = m_simulationShader->addCacheableShaderFromSourceFile(QOpenGLShader::ShaderTypeBit::Compute, ":/shaders/raytrace.glsl");
     if (raytraceShaderReadSucceeded == false)
     {
@@ -294,12 +336,30 @@ void SimulationEngine::initializeShaders()
         throw std::runtime_error(m_skyShader->log().toUtf8());
     }
     qInfo("Sky shader program compilation and linking successful");
+
+    qInfo("Initializing guide marking shader");
+    m_guideShader = new QOpenGLShaderProgram(this);
+    bool guideShaderReadSucceeded = m_guideShader->addCacheableShaderFromSourceFile(QOpenGLShader::ShaderTypeBit::Compute, ":/shaders/guide.glsl");
+    if (guideShaderReadSucceeded == false)
+    {
+        qWarning("Reading guide marking shader failed");
+        throw std::runtime_error(m_guideShader->log().toUtf8());
+    }
+    qInfo("Guide marking shader successfully initialized");
+
+    if (m_guideShader->link() == false)
+    {
+        qWarning("Compiling and linking guide marking shader failed");
+        throw std::runtime_error(m_guideShader->log().toUtf8());
+    }
+    qInfo("Guide marking shader program compilation and linking successful");
 }
 
 void SimulationEngine::initializeTextures()
 {
     m_simulationTexture = std::make_unique<OpenGL::Texture>(m_outputWidth, m_outputHeight, 0, OpenGL::TextureType::Color);
-    m_backgroundTexture = std::make_unique<OpenGL::Texture>(m_outputWidth, m_outputHeight, 2, OpenGL::TextureType::Color);
+    m_backgroundTexture = std::make_unique<OpenGL::Texture>(m_outputWidth, m_outputHeight, 1, OpenGL::TextureType::Color);
+    m_guideTexture = std::make_unique<OpenGL::Texture>(m_outputWidth, m_outputHeight, 2, OpenGL::TextureType::Monochrome);
 }
 
 void SimulationEngine::resizeOutputTextureCallback(const unsigned int width, const unsigned int height)
@@ -309,6 +369,7 @@ void SimulationEngine::resizeOutputTextureCallback(const unsigned int width, con
 
     m_simulationTexture.reset();
     m_backgroundTexture.reset();
+    m_guideTexture.reset();
 
     initializeTextures();
     clear();
diff --git a/src/haloray-core/simulation/simulationEngine.h b/src/haloray-core/simulation/simulationEngine.h
index 62ba29d..e7b739b 100644
--- a/src/haloray-core/simulation/simulationEngine.h
+++ b/src/haloray-core/simulation/simulationEngine.h
@@ -41,6 +41,9 @@ class SimulationEngine : public QObject, protected QOpenGLFunctions_4_4_Core
     Atmosphere getAtmosphere() const;
     void setAtmosphere(Atmosphere);
 
+    bool getGuidesEnabled() const;
+    void setGuidesEnabled(bool);
+
     void lockCameraToLightSource(bool locked);
 
     void setMultipleScatteringProbability(double);
@@ -48,6 +51,7 @@ class SimulationEngine : public QObject, protected QOpenGLFunctions_4_4_Core
 
     unsigned int getOutputTextureHandle() const;
     unsigned int getBackgroundTextureHandle() const;
+    unsigned int getGuideTextureHandle() const;
 
     void resizeOutputTextureCallback(const unsigned int width, const unsigned int height);
 
@@ -58,6 +62,7 @@ class SimulationEngine : public QObject, protected QOpenGLFunctions_4_4_Core
     void atmosphereChanged(Atmosphere);
     void lockCameraToLightSourceChanged(bool);
     void multipleScatteringProbabilityChanged(double);
+    void guidesToggled(bool);
 
 private:
     void initializeShaders();
@@ -68,10 +73,12 @@ class SimulationEngine : public QObject, protected QOpenGLFunctions_4_4_Core
     unsigned int m_outputHeight;
     std::mt19937 m_mersenneTwister;
     std::uniform_int_distribution<unsigned int> m_uniformDistribution;
-    std::unique_ptr<QOpenGLShaderProgram> m_simulationShader;
+    QOpenGLShaderProgram *m_simulationShader;
     QOpenGLShaderProgram *m_skyShader;
+    QOpenGLShaderProgram *m_guideShader;
     std::unique_ptr<OpenGL::Texture> m_simulationTexture;
     std::unique_ptr<OpenGL::Texture> m_backgroundTexture;
+    std::unique_ptr<OpenGL::Texture> m_guideTexture;
 
     Camera m_camera;
     LightSource m_light;
@@ -85,6 +92,7 @@ class SimulationEngine : public QObject, protected QOpenGLFunctions_4_4_Core
     std::shared_ptr<CrystalPopulationRepository> m_crystalRepository;
     float m_sunSpectrumCache[31];
     Atmosphere m_atmosphere;
+    bool m_guidesEnabled;
 };
 
 }