<p>There are a number of minor problems with the questioner's (the OP, original poster after this), but they aren't the main source of the problem. The biggest problem is the slant in the plot. I'll go after the lesser problems first.</p>
<p>From a comment I made,</p>
<blockquote>
<p>One thing you are definitely doing wrong is mixing and matching J2000 and ECEF. What you are labeling as latitude and longitude should be declination and right ascension. To use latitude and longitude, you should be transforming <code>observerPos</code> and <code>zenith</code> to J2000. You'll need a transformation matrix to do that. SPICE offers that capability via <code>pxform_c</code>. The transformation matrix from ITRF93 to J2000 does just that.</p>
</blockquote>
<p>Do not mix and match reference frames in SPICE (or in any other technical software package capable of handling multiple frames of reference).</p>
<p>A lesser problem is using <code>irradiance = SOLAR_CONSTANT * sin(altitude);</code>. Since the code in question is attempting to calculate solar intensity at sea level, a better approach is to calculate air mass, and from that, calculate intensity.</p>
<p>Neither of these address the slant issue. Solar noon (in places of the world where they have a daily solar noon) occurs at the same time regardless of latitude (but definitely dependent on longitude). Latitude affects how high the Sun will be over the horizon (and whether it is over the horizon). The image should not have a slant.</p>
<p>The problem turned out to be using <code>std::vector&lt;std::vector&lt;std::vector&lt;double&gt;&gt;&gt;</code> for <code>solarIrradiance</code> and <code>luminance</code> with an initial value that populated those matrices. For example,</p>
<pre><code>std::vector&lt;std::vector&lt;std::vector&lt;double&gt;&gt;&gt;
    solarIrradiance(times.size(),
        std::vector&lt;std::vector&lt;double&gt;&gt;(latitudes.size(),
            std::vector&lt;double&gt;(longitudes.size(), 0.0)));
</code></pre>
<p>And similar for <code>luminance</code>. The C++ standard does not guarantee that this will be constructed as a contiguous 3D row major array. In particular, passing <code>&amp;solarIrradiance[0][0][0]</code> (and similar for <code>luminance</code>) to a NetCDF function that expects the data pointer to point to the first element of a contiguous row major order multidimensional array might not work. This is yet another <em>wonderful</em> example of undefined behavior in C++ (<em>wonderful</em> written rather sarcastically). In general, when I worked (I am now retired), I strongly discouraged using <code>std::vector</code> for fixed-size multidimensional arrays. I also discouraged using <code>std::array</code>, but for a different reason: multidimensional <code>std::array</code> is column major rather than row major.</p>
<p>What I did encourage when I did work, and what the OP eventually did do, was to write a three dimensional matrix class (or class template) that overloads the function call operator (<code>Class::operator()</code>) that takes three integer type arguments that specify the desired element in row major order. Two versions of the overload are often useful, a non-const version that returns a reference to the desired element, and a const version that returns the value of that desired element. This worked. No more slant!</p>
<p><a href="https://i.sstatic.net/FQga7NVo.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/FQga7NVo.png" alt="Worldwide solar luminance at the start time of the time period of interest" /></a><br />
<sub>Source: Image created by NitroxHead, the author of the question, after changing the use of a vector of a vector of a vector to a contiguous 3D row major array</sub></p>
<p>There's a lesson to be learned here: Don't invoke undefined behavior. That is easier said than done because there is so much undefined behavior in C++ that the standard doesn't bother to list them. (The C standard does, or at least tries to. The last iteration of the C standard I read had a 14 page long appendix listing undefined behavior in C. The list was acknowledged to most likely be incomplete. But at least they tried. The C++ standard doesn't even try.)</p>

CC BY-SA 4.0
2024-09-29T05:45:42.967Z
2024-09-29T08:26:07.157Z
2024-09-29T08:26:07.157Z
David Hammen
<p>I'm a noob at astronomy, but I needed to calculate a very simple theoretical solar irradiance and luminance. With the help of LLMs, I came up with a C++ script, but the output doesn't look correct. I'm open to any input specifically seeking any help to fix it!</p>
<p>Script calculates solar irradiance and luminance based on time, latitude, and longitude using the SPICE toolkit for celestial mechanics, and outputs the data to a NetCDF file. My code compiles and runs without errors, but the output for solar irradiance and luminance seems incorrect. I suspect the problem lies either in the calculation logic or how the data is stored.</p>
<p>Issue:</p>
<p>I don't know how exactly I can but it in the words but can be seen in the image.</p>
<p>My theory is as different locations are calculated time shifts but I'm no expert.</p>
<p><a href="https://i.sstatic.net/JpYT9UF2.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JpYT9UF2.png" alt="output values" /></a></p>
<p>Here is the code;</p>
<pre><code>#include &lt;iostream&gt;
#include &lt;vector&gt;
#include &lt;cmath&gt;
#include &lt;string&gt;
#include &lt;chrono&gt;
#include &lt;netcdf&gt;
#include &quot;./include/SpiceUsr.h&quot;

// Constants
const double SOLAR_CONSTANT = 1361.0; // W/m²
const double CONVERSION_FACTOR = 93.0; // lumens/m² per W/m²
const double DEG_TO_RAD = M_PI / 180.0;

void calculateSolarData(double et, double lat, double lon, double&amp; irradiance, double&amp; luminance) {
    double sunPos[3], sunLt;
    
    // Get the Sun's position relative to Earth
    spkpos_c(&quot;SUN&quot;, et, &quot;J2000&quot;, &quot;LT+S&quot;, &quot;EARTH&quot;, sunPos, &amp;sunLt);

    // Convert lat/lon to rectangular coordinates on Earth's surface
    double observerPos[3];
    double alt = 0.0;  // Assume observer is at sea level
    georec_c(lon * DEG_TO_RAD, lat * DEG_TO_RAD, alt, 6378.137, 0.0033528131, observerPos);

    // Calculate the vector from observer to Sun
    double observerToSun[3];
    vsub_c(sunPos, observerPos, observerToSun);

    // Calculate the local zenith vector
    double zenith[3];
    surfnm_c(6378.137, 6378.137, 6356.752, observerPos, zenith);

    // Calculate the angle between zenith and Sun direction
    double cosSolarZenithAngle = vdot_c(zenith, observerToSun) / vnorm_c(observerToSun);
    double solarZenithAngle = acos(cosSolarZenithAngle);
    double altitude = halfpi_c() - solarZenithAngle;

    if (altitude &gt; 0) {
        irradiance = SOLAR_CONSTANT * sin(altitude);
        luminance = irradiance * CONVERSION_FACTOR;
    } else {
        irradiance = 0.0;
        luminance = 0.0;
    }
}

// Function to create NetCDF output
void createNetCDFOutput(const std::string&amp; filename,
                        const std::vector&lt;double&gt;&amp; times,
                        const std::vector&lt;double&gt;&amp; latitudes,
                        const std::vector&lt;double&gt;&amp; longitudes,
                        const std::vector&lt;std::vector&lt;std::vector&lt;double&gt;&gt;&gt;&amp; solarIrradiance,
                        const std::vector&lt;std::vector&lt;std::vector&lt;double&gt;&gt;&gt;&amp; luminance) {
    try {
        // Create NetCDF file
        netCDF::NcFile dataFile(filename, netCDF::NcFile::replace);

        // Define dimensions
        auto timeDim = dataFile.addDim(&quot;time&quot;, times.size());
        auto latDim = dataFile.addDim(&quot;latitude&quot;, latitudes.size());
        auto lonDim = dataFile.addDim(&quot;longitude&quot;, longitudes.size());

        // Define coordinate variables
        auto timeVar = dataFile.addVar(&quot;time&quot;, netCDF::ncDouble, timeDim);
        auto latVar = dataFile.addVar(&quot;latitude&quot;, netCDF::ncDouble, latDim);
        auto lonVar = dataFile.addVar(&quot;longitude&quot;, netCDF::ncDouble, lonDim);

        // Define data variables
        std::vector&lt;netCDF::NcDim&gt; dims = {timeDim, latDim, lonDim};
        auto irradianceVar = dataFile.addVar(&quot;solar_irradiance&quot;, netCDF::ncDouble, dims);
        auto luminanceVar = dataFile.addVar(&quot;luminance&quot;, netCDF::ncDouble, dims);

        // Add attributes
        timeVar.putAtt(&quot;units&quot;, &quot;hours since 1950-01-01 00:00:00&quot;);
        timeVar.putAtt(&quot;calendar&quot;, &quot;gregorian&quot;);

        irradianceVar.putAtt(&quot;units&quot;, &quot;W/m²&quot;);
        irradianceVar.putAtt(&quot;long_name&quot;, &quot;Solar Irradiance&quot;);

        luminanceVar.putAtt(&quot;units&quot;, &quot;lux&quot;);
        luminanceVar.putAtt(&quot;long_name&quot;, &quot;Luminance&quot;);

        // Global attributes
        auto nowTime = std::chrono::system_clock::now();
        auto now = std::chrono::system_clock::to_time_t(nowTime);
        dataFile.putAtt(&quot;Conventions&quot;, &quot;CF-1.6&quot;);
        dataFile.putAtt(&quot;title&quot;, &quot;Solar Irradiance and Luminance Data&quot;);
        dataFile.putAtt(&quot;source&quot;, &quot;Generated from SPICE calculations&quot;);
        dataFile.putAtt(&quot;constants&quot;, &quot;93 # lumens/m² per W/m² \n 1361 # W/m²&quot;);
        dataFile.putAtt(&quot;MadeBy&quot;, &quot;NitroxHead&quot;);
        dataFile.putAtt(&quot;history&quot;, &quot;Created on &quot; + std::string(std::ctime(&amp;now)));

        // Write data
        timeVar.putVar(times.data());
        latVar.putVar(latitudes.data());
        lonVar.putVar(longitudes.data());

        // Write irradiance and luminance data
        std::vector&lt;size_t&gt; startp = {0, 0, 0};
        std::vector&lt;size_t&gt; countp = {times.size(), latitudes.size(), longitudes.size()};
        irradianceVar.putVar(startp, countp, &amp;solarIrradiance[0][0][0]);
        luminanceVar.putVar(startp, countp, &amp;luminance[0][0][0]);

    } catch (netCDF::exceptions::NcException&amp; e) {
        std::cerr &lt;&lt; &quot;NetCDF exception: &quot; &lt;&lt; e.what() &lt;&lt; std::endl;
        exit(1);
    }

    std::cout &lt;&lt; &quot;Created &quot; &lt;&lt; filename &lt;&lt; std::endl;
}

int main() {
    // Load SPICE kernels
    furnsh_c(&quot;de421.bsp&quot;);
    furnsh_c(&quot;naif0012.tls&quot;);
    furnsh_c(&quot;pck00010.tpc&quot;);

    // Define the datetime range with more precise format
    std::string startDate = &quot;2019-01-01T00:00:00.000&quot;;
    std::string endDate = &quot;2019-01-02T00:00:00.000&quot;;

    double startEt, endEt;
    str2et_c(startDate.c_str(), &amp;startEt);
    str2et_c(endDate.c_str(), &amp;endEt);

    // Print loaded start and end times in ET for debug
    std::cout &lt;&lt; &quot;Start ET: &quot; &lt;&lt; startEt &lt;&lt; &quot;, End ET: &quot; &lt;&lt; endEt &lt;&lt; std::endl;

    // Convert ET to UTC and print it for verification
    char utcStr[30];
    et2utc_c(startEt, &quot;ISOC&quot;, 3, 30, utcStr);
    std::cout &lt;&lt; &quot;Start Time (UTC): &quot; &lt;&lt; utcStr &lt;&lt; std::endl;
    et2utc_c(endEt, &quot;ISOC&quot;, 3, 30, utcStr);
    std::cout &lt;&lt; &quot;End Time (UTC): &quot; &lt;&lt; utcStr &lt;&lt; std::endl;

    // Create a grid of latitude and longitude coordinates
    std::vector&lt;double&gt; latitudes, longitudes;
    for (double lat = -90.0; lat &lt; 90.0; lat += 0.08333588) {
        latitudes.push_back(lat);
    }
    for (double lon = -180.0; lon &lt; 180.0; lon += 0.08333588) {
        longitudes.push_back(lon);
    }

    // Generate time steps (8 intervals of 3 hours)
    std::vector&lt;double&gt; times;
    double timeStep = 3 * 3600; // 3 hours in seconds
    for (int i = 0; i &lt; 8; ++i) {
        times.push_back(startEt + i * timeStep);
    }

    // Print generated times in UTC for debug
    for (double t : times) {
        et2utc_c(t, &quot;C&quot;, 3, 30, utcStr);
        std::cout &lt;&lt; &quot;Time (UTC): &quot; &lt;&lt; utcStr &lt;&lt; std::endl;
    }

    // Initialize data structures
    std::vector&lt;std::vector&lt;std::vector&lt;double&gt;&gt;&gt; solarIrradiance(times.size(),
        std::vector&lt;std::vector&lt;double&gt;&gt;(latitudes.size(),
        std::vector&lt;double&gt;(longitudes.size(), 0.0)));
    std::vector&lt;std::vector&lt;std::vector&lt;double&gt;&gt;&gt; luminance(times.size(),
        std::vector&lt;std::vector&lt;double&gt;&gt;(latitudes.size(),
        std::vector&lt;double&gt;(longitudes.size(), 0.0)));

    // Calculate solar data
    for (size_t t = 0; t &lt; times.size(); ++t) {
        for (size_t i = 0; i &lt; latitudes.size(); ++i) {
            for (size_t j = 0; j &lt; longitudes.size(); ++j) {
                double irr, lum;
                calculateSolarData(times[t], latitudes[i], longitudes[j], irr, lum);
                solarIrradiance[t][i][j] = irr;
                luminance[t][i][j] = lum;
            }
        }
    }

    // Create NetCDF output file
    createNetCDFOutput(&quot;solar_data_20190101.nc&quot;, times, latitudes, longitudes, solarIrradiance, luminance);

    // Unload SPICE kernels
    kclear_c();

    return 0;
}
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2024-09-23T16:05:54.567Z
What kind of output were you expecting?  What you have seems plausible at first glance.
Greg Miller
CC BY-SA 4.0
2024-09-24T07:14:03.527Z
@GregMiller I got some result from another script and result looked like this: [link](https://puntocopy.com/nitroxhead/327d182c366c497882f69236db3a7262/lumi2_py_prev1.png). I believe it supposed to look like the last example.
NitroxHead
CC BY-SA 4.0
2024-09-24T14:03:59.410Z
This definitely does not look right. That strong tilt should not be there. Figure out what is causing that will fix a lot of problems. I can't see what's wrong. One thing you are definitely doing wrong is mixing and matching J2000 and ECEF. What you are labeling as latitude and longitude should be declination and right ascension. To use latitude and longitude, you should be transforming `observerPos` and `zenith` to J2000. You'll need a transformation matrix to do that. SPICE offers that capability via `pxform_c`. The transformation matrix from ITRF93 to J2000 does just that.
David Hammen
CC BY-SA 4.0
2024-09-24T14:10:22.440Z
Another thing you are doing wrong is `irradiance = SOLAR_CONSTANT * sin(altitude);` Did the LLM give you that expression? That's just wrong. A better model would be to calculate air mass at sea level and from that calculate the intensity of sunlight.
David Hammen
CC BY-SA 4.0
2024-09-24T15:16:08.647Z
@DavidHammen Thanks for the suggestions. I corrected the function and looks like there is still something wrong (but now more precisely :) ). [current function](https://writexo.com/2gdfm7tc) [current image](https://ibb.co/cgKZ4Xx).

As for calculation of `irradiance`, yes. I know it is quite wrong, but there is a bug somewhere before that, want to fix that before applying more complex model.
NitroxHead
CC BY-SA 4.0
2024-09-24T15:34:56.440Z
@NitroxHead That tilt is still present. It shouldn't be there. Look at time zones. The boundaries more or less follow lines of longitude. (They don't because of political concerns). There is no tilt, however. I can't see what is causing that tilt.
David Hammen
CC BY-SA 4.0
2024-09-24T15:56:38.517Z
But at least that tilted blob should now be marching across the globe over your one day time span.
David Hammen
CC BY-SA 4.0
2024-09-24T19:32:19.977Z
I suspect the root of the problem might be with the use of `std::vector<std::vector<std::vector<double>>>` for `solarIrradiance` and `luminance`. Using `std::vector` for a fixed size multidimensional array is in general not a good idea. While each of the innermost vectors will be contiguous, there is no guarantee of contiguousness regarding the outer elements. There is no guarantee that (for example) `&luminance[0][0][0]` will access data in the order you expect. (continued)
David Hammen
CC BY-SA 4.0
2024-09-24T19:41:36.797Z
I suggest writing a class whose constructor takes the three dimensions as arguments, stores those as data members (for short, `l`. `m`, and `n`) and creates a `std::vector` of size `l*m*n`. Override the function operator to serve as the multidimensional index into the vector that calculates the index into the vector in row-major order, and finally, provides a function `data` that gives access to the underlying 1D vector as a `double*` (or `const double*`). Note that the underlying data is guaranteed to be contiguous because that's how `std::vector` works, and is in the proper order.
David Hammen
CC BY-SA 4.0
2024-09-24T19:47:07.450Z
For example, `double& 3DMatrix::operator () (int i, int j, int k) { return data_[i*m*n + j*n + k]; }`. You might want a `const operator()` as well that returns a `double`. Note that there is no protection against index out of bounds. Then again, the same applies to `std::vector`.
David Hammen
CC BY-SA 4.0
2024-09-29T01:54:45.647Z
@DavidHammen I applied your suggestion. I couldn't believe my eyes it worked!
NitroxHead
CC BY-SA 4.0
2024-09-29T02:02:49.480Z
@DavidHammen Here is what it looks like [image](https://ibb.co/t3pr7Hs) .I would have continued wasting my hours looking for mistakes in the calculation function if it weren't for you. Thank you! This is my first time using Stack Exchange I couldn't figure out how to accept your comment and answer... I will share the code on GitHub addressing your contribution on commit description.
NitroxHead
CC BY-SA 4.0
2024-09-29T05:46:49.290Z
@NitroxHead You couldn't accept my answer because comments are not answers. I have now written up a summary of this long comment section as an answer.
David Hammen
<the-sun><data-analysis><physics-simulations>
Incorrect Solar Data Output Using SPICE Toolkit
NitroxHead
# Incorrect Solar Data Output Using SPICE Toolkit

(Asked by: username_0 on 2024-09-23)

I'm a noob at astronomy, but I needed to calculate a very simple theoretical solar irradiance and luminance. With the help of LLMs, I came up with a C++ script, but the output doesn't look correct. I'm open to any input specifically seeking any help to fix it!

Script calculates solar irradiance and luminance based on time, latitude, and longitude using the SPICE toolkit for celestial mechanics, and outputs the data to a NetCDF file. My code compiles and runs without errors, but the output for solar irradiance and luminance seems incorrect. I suspect the problem lies either in the calculation logic or how the data is stored.

Issue:

I don't know how exactly I can but it in the words but can be seen in the image.

My theory is as different locations are calculated time shifts but I'm no expert.

[![output values](https://i.sstatic.net/JpYT9UF2.png)](https://i.sstatic.net/JpYT9UF2.png)

Here is the code;

```
#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include <chrono>
#include <netcdf>
#include "./include/SpiceUsr.h"

// Constants
const double SOLAR_CONSTANT = 1361.0; // W/m²
const double CONVERSION_FACTOR = 93.0; // lumens/m² per W/m²
const double DEG_TO_RAD = M_PI / 180.0;

void calculateSolarData(double et, double lat, double lon, double& irradiance, double& luminance) {
    double sunPos[3], sunLt;
    
    // Get the Sun's position relative to Earth
    spkpos_c("SUN", et, "J2000", "LT+S", "EARTH", sunPos, &sunLt);

    // Convert lat/lon to rectangular coordinates on Earth's surface
    double observerPos[3];
    double alt = 0.0;  // Assume observer is at sea level
    georec_c(lon * DEG_TO_RAD, lat * DEG_TO_RAD, alt, 6378.137, 0.0033528131, observerPos);

    // Calculate the vector from observer to Sun
    double observerToSun[3];
    vsub_c(sunPos, observerPos, observerToSun);

    // Calculate the local zenith vector
    double zenith[3];
    surfnm_c(6378.137, 6378.137, 6356.752, observerPos, zenith);

    // Calculate the angle between zenith and Sun direction
    double cosSolarZenithAngle = vdot_c(zenith, observerToSun) / vnorm_c(observerToSun);
    double solarZenithAngle = acos(cosSolarZenithAngle);
    double altitude = halfpi_c() - solarZenithAngle;

    if (altitude > 0) {
        irradiance = SOLAR_CONSTANT * sin(altitude);
        luminance = irradiance * CONVERSION_FACTOR;
    } else {
        irradiance = 0.0;
        luminance = 0.0;
    }
}

// Function to create NetCDF output
void createNetCDFOutput(const std::string& filename,
                        const std::vector<double>& times,
                        const std::vector<double>& latitudes,
                        const std::vector<double>& longitudes,
                        const std::vector<std::vector<std::vector<double>>>& solarIrradiance,
                        const std::vector<std::vector<std::vector<double>>>& luminance) {
    try {
        // Create NetCDF file
        netCDF::NcFile dataFile(filename, netCDF::NcFile::replace);

        // Define dimensions
        auto timeDim = dataFile.addDim("time", times.size());
        auto latDim = dataFile.addDim("latitude", latitudes.size());
        auto lonDim = dataFile.addDim("longitude", longitudes.size());

        // Define coordinate variables
        auto timeVar = dataFile.addVar("time", netCDF::ncDouble, timeDim);
        auto latVar = dataFile.addVar("latitude", netCDF::ncDouble, latDim);
        auto lonVar = dataFile.addVar("longitude", netCDF::ncDouble, lonDim);

        // Define data variables
        std::vector<netCDF::NcDim> dims = {timeDim, latDim, lonDim};
        auto irradianceVar = dataFile.addVar("solar_irradiance", netCDF::ncDouble, dims);
        auto luminanceVar = dataFile.addVar("luminance", netCDF::ncDouble, dims);

        // Add attributes
        timeVar.putAtt("units", "hours since 1950-01-01 00:00:00");
        timeVar.putAtt("calendar", "gregorian");

        irradianceVar.putAtt("units", "W/m²");
        irradianceVar.putAtt("long_name", "Solar Irradiance");

        luminanceVar.putAtt("units", "lux");
        luminanceVar.putAtt("long_name", "Luminance");

        // Global attributes
        auto nowTime = std::chrono::system_clock::now();
        auto now = std::chrono::system_clock::to_time_t(nowTime);
        dataFile.putAtt("Conventions", "CF-1.6");
        dataFile.putAtt("title", "Solar Irradiance and Luminance Data");
        dataFile.putAtt("source", "Generated from SPICE calculations");
        dataFile.putAtt("constants", "93 # lumens/m² per W/m² \n 1361 # W/m²");
        dataFile.putAtt("MadeBy", "username_0");
        dataFile.putAtt("history", "Created on " + std::string(std::ctime(&now)));

        // Write data
        timeVar.putVar(times.data());
        latVar.putVar(latitudes.data());
        lonVar.putVar(longitudes.data());

        // Write irradiance and luminance data
        std::vector<size_t> startp = {0, 0, 0};
        std::vector<size_t> countp = {times.size(), latitudes.size(), longitudes.size()};
        irradianceVar.putVar(startp, countp, &solarIrradiance[0][0][0]);
        luminanceVar.putVar(startp, countp, &luminance[0][0][0]);

    } catch (netCDF::exceptions::NcException& e) {
        std::cerr << "NetCDF exception: " << e.what() << std::endl;
        exit(1);
    }

    std::cout << "Created " << filename << std::endl;
}

int main() {
    // Load SPICE kernels
    furnsh_c("de421.bsp");
    furnsh_c("naif0012.tls");
    furnsh_c("pck00010.tpc");

    // Define the datetime range with more precise format
    std::string startDate = "2019-01-01T00:00:00.000";
    std::string endDate = "2019-01-02T00:00:00.000";

    double startEt, endEt;
    str2et_c(startDate.c_str(), &startEt);
    str2et_c(endDate.c_str(), &endEt);

    // Print loaded start and end times in ET for debug
    std::cout << "Start ET: " << startEt << ", End ET: " << endEt << std::endl;

    // Convert ET to UTC and print it for verification
    char utcStr[30];
    et2utc_c(startEt, "ISOC", 3, 30, utcStr);
    std::cout << "Start Time (UTC): " << utcStr << std::endl;
    et2utc_c(endEt, "ISOC", 3, 30, utcStr);
    std::cout << "End Time (UTC): " << utcStr << std::endl;

    // Create a grid of latitude and longitude coordinates
    std::vector<double> latitudes, longitudes;
    for (double lat = -90.0; lat < 90.0; lat += 0.08333588) {
        latitudes.push_back(lat);
    }
    for (double lon = -180.0; lon < 180.0; lon += 0.08333588) {
        longitudes.push_back(lon);
    }

    // Generate time steps (8 intervals of 3 hours)
    std::vector<double> times;
    double timeStep = 3 * 3600; // 3 hours in seconds
    for (int i = 0; i < 8; ++i) {
        times.push_back(startEt + i * timeStep);
    }

    // Print generated times in UTC for debug
    for (double t : times) {
        et2utc_c(t, "C", 3, 30, utcStr);
        std::cout << "Time (UTC): " << utcStr << std::endl;
    }

    // Initialize data structures
    std::vector<std::vector<std::vector<double>>> solarIrradiance(times.size(),
        std::vector<std::vector<double>>(latitudes.size(),
        std::vector<double>(longitudes.size(), 0.0)));
    std::vector<std::vector<std::vector<double>>> luminance(times.size(),
        std::vector<std::vector<double>>(latitudes.size(),
        std::vector<double>(longitudes.size(), 0.0)));

    // Calculate solar data
    for (size_t t = 0; t < times.size(); ++t) {
        for (size_t i = 0; i < latitudes.size(); ++i) {
            for (size_t j = 0; j < longitudes.size(); ++j) {
                double irr, lum;
                calculateSolarData(times[t], latitudes[i], longitudes[j], irr, lum);
                solarIrradiance[t][i][j] = irr;
                luminance[t][i][j] = lum;
            }
        }
    }

    // Create NetCDF output file
    createNetCDFOutput("solar_data_20190101.nc", times, latitudes, longitudes, solarIrradiance, luminance);

    // Unload SPICE kernels
    kclear_c();

    return 0;
}
```

```

**Tags:** the-sun, data-analysis, physics-simulations

> **Comments:**
> **username_1**: What kind of output were you expecting?  What you have seems plausible at first glance.  
> **username_0**: @username_1 I got some result from another script and result looked like this: [link](https://puntocopy.com/username_0/327d182c366c497882f69236db3a7262/lumi2_py_prev1.png). I believe it supposed to look like the last example.  
> **username_2**: This definitely does not look right. That strong tilt should not be there. Figure out what is causing that will fix a lot of problems. I can't see what's wrong. One thing you are definitely doing wrong is mixing and matching J2000 and ECEF. What you are labeling as latitude and longitude should be declination and right ascension. To use latitude and longitude, you should be transforming `observerPos` and `zenith` to J2000. You'll need a transformation matrix to do that. SPICE offers that capability via `pxform_c`. The transformation matrix from ITRF93 to J2000 does just that.  
> **username_2**: Another thing you are doing wrong is `irradiance = SOLAR_CONSTANT * sin(altitude);` Did the LLM give you that expression? That's just wrong. A better model would be to calculate air mass at sea level and from that calculate the intensity of sunlight.  
> **username_0**: @username_2 Thanks for the suggestions. I corrected the function and looks like there is still something wrong (but now more precisely :) ). [current function](https://writexo.com/2gdfm7tc) [current image](https://ibb.co/cgKZ4Xx).  As for calculation of `irradiance`, yes. I know it is quite wrong, but there is a bug somewhere before that, want to fix that before applying more complex model.  
> **username_2**: @username_0 That tilt is still present. It shouldn't be there. Look at time zones. The boundaries more or less follow lines of longitude. (They don't because of political concerns). There is no tilt, however. I can't see what is causing that tilt.  
> **username_2**: But at least that tilted blob should now be marching across the globe over your one day time span.  
> **username_2**: I suspect the root of the problem might be with the use of `std::vector<std::vector<std::vector<double>>>` for `solarIrradiance` and `luminance`. Using `std::vector` for a fixed size multidimensional array is in general not a good idea. While each of the innermost vectors will be contiguous, there is no guarantee of contiguousness regarding the outer elements. There is no guarantee that (for example) `&luminance[0][0][0]` will access data in the order you expect. (continued)  
> **username_2**: I suggest writing a class whose constructor takes the three dimensions as arguments, stores those as data members (for short, `l`. `m`, and `n`) and creates a `std::vector` of size `l*m*n`. Override the function operator to serve as the multidimensional index into the vector that calculates the index into the vector in row-major order, and finally, provides a function `data` that gives access to the underlying 1D vector as a `double*` (or `const double*`). Note that the underlying data is guaranteed to be contiguous because that's how `std::vector` works, and is in the proper order.  
> **username_2**: For example, `double& 3DMatrix::operator () (int i, int j, int k) { return data_[i*m*n + j*n + k]; }`. You might want a `const operator()` as well that returns a `double`. Note that there is no protection against index out of bounds. Then again, the same applies to `std::vector`.  
> **username_0**: @username_2 I applied your suggestion. I couldn't believe my eyes it worked!  
> **username_0**: @username_2 Here is what it looks like [image](https://ibb.co/t3pr7Hs) .I would have continued wasting my hours looking for mistakes in the calculation function if it weren't for you. Thank you! This is my first time using Stack Exchange I couldn't figure out how to accept your comment and answer... I will share the code on GitHub addressing your contribution on commit description.  
> **username_2**: @username_0 You couldn't accept my answer because comments are not answers. I have now written up a summary of this long comment section as an answer.  

---

## 1 answers

### Answer by username_2 on 2024-09-29. Score: 2 | Accepted answer

There are a number of minor problems with the questioner's (the OP, original poster after this), but they aren't the main source of the problem. The biggest problem is the slant in the plot. I'll go after the lesser problems first.

From a comment I made,

> One thing you are definitely doing wrong is mixing and matching J2000 and ECEF. What you are labeling as latitude and longitude should be declination and right ascension. To use latitude and longitude, you should be transforming `observerPos` and `zenith` to J2000. You'll need a transformation matrix to do that. SPICE offers that capability via `pxform_c`. The transformation matrix from ITRF93 to J2000 does just that.

Do not mix and match reference frames in SPICE (or in any other technical software package capable of handling multiple frames of reference).

A lesser problem is using `irradiance = SOLAR_CONSTANT * sin(altitude);`. Since the code in question is attempting to calculate solar intensity at sea level, a better approach is to calculate air mass, and from that, calculate intensity.

Neither of these address the slant issue. Solar noon (in places of the world where they have a daily solar noon) occurs at the same time regardless of latitude (but definitely dependent on longitude). Latitude affects how high the Sun will be over the horizon (and whether it is over the horizon). The image should not have a slant.

The problem turned out to be using `std::vector<std::vector<std::vector<double>>>` for `solarIrradiance` and `luminance` with an initial value that populated those matrices. For example,

```
std::vector<std::vector<std::vector<double>>>
    solarIrradiance(times.size(),
        std::vector<std::vector<double>>(latitudes.size(),
            std::vector<double>(longitudes.size(), 0.0)));

```

And similar for `luminance`. The C++ standard does not guarantee that this will be constructed as a contiguous 3D row major array. In particular, passing `&solarIrradiance[0][0][0]` (and similar for `luminance`) to a NetCDF function that expects the data pointer to point to the first element of a contiguous row major order multidimensional array might not work. This is yet another *wonderful* example of undefined behavior in C++ (*wonderful* written rather sarcastically). In general, when I worked (I am now retired), I strongly discouraged using `std::vector` for fixed-size multidimensional arrays. I also discouraged using `std::array`, but for a different reason: multidimensional `std::array` is column major rather than row major.

What I did encourage when I did work, and what the OP eventually did do, was to write a three dimensional matrix class (or class template) that overloads the function call operator (`Class::operator()`) that takes three integer type arguments that specify the desired element in row major order. Two versions of the overload are often useful, a non-const version that returns a reference to the desired element, and a const version that returns the value of that desired element. This worked. No more slant!

[![Worldwide solar luminance at the start time of the time period of interest](https://i.sstatic.net/FQga7NVo.png)](https://i.sstatic.net/FQga7NVo.png)  
Source: Image created by username_0, the author of the question, after changing the use of a vector of a vector of a vector to a contiguous 3D row major array

There's a lesson to be learned here: Don't invoke undefined behavior. That is easier said than done because there is so much undefined behavior in C++ that the standard doesn't bother to list them. (The C standard does, or at least tries to. The last iteration of the C standard I read had a 14 page long appendix listing undefined behavior in C. The list was acknowledged to most likely be incomplete. But at least they tried. The C++ standard doesn't even try.)


---

