close-circle
Close
0%
0%

Daylight Geolocator Remix

This is a remix of https://hackaday.io/project/28550-light-level-geolocator

Similar projects worth following
close
A remix of https://hackaday.io/project/28550-light-level-geolocator.A device that logs daylight and tries to workout your latitude and longitude.Might have been useful before the days of GPS.

Project Stages

  • The mathematics:
    • Calculating modified Julian Day number (JD)
    • Calculating Greenwich Mean Sidereal Time (GMST)
    • Calculating Universal Time
    • Calculating the Right Ascension (RA) and Declination (Dec) of the Sun.
    • Calculating sun rise and sun set or any other position of the sun.
    • Calculating your Longitude based on sun rise and sun set times.
    • Calculating your Latitude based on sun rise and sun set times.
    • Some code to try.
  • A box to put the electronic in.
  • What could be so hard with samping daylight:
    • The need for an accurate Real Time Clock
    • The rotary encoder as an input device
    • Code to read the rotary encoder
    • A menu system
    • Setting the time
    • Saving your data
    • Remembering settings
    • Downloading your data
    • Dealing with power outages (restarting as if nothing happened)
  •  Processing your data:
    • TBC ...

JPEG Image - 109.21 kB - 02/04/2018 at 11:46

eye
Preview
download-circle
Download

Adobe Portable Document Format - 85.09 kB - 02/04/2018 at 11:46

eye
Preview
download-circle
Download

JPEG Image - 35.79 kB - 02/04/2018 at 11:46

eye
Preview
download-circle
Download

Portable Network Graphics (PNG) - 6.90 kB - 02/04/2018 at 11:46

eye
Preview
download-circle
Download

x-wine-extension-bas - 30.63 kB - 02/04/2018 at 11:46

download-circle
Download

View all 26 files

  • PROJECT REBOOT

    agp.cooper5 days ago 0 comments

    Project Reboot - A major Change Direction

    I have given up on accurate estimation of sunrise (or dawn) and sunset (or dusk) directly from daylight measurements. It seems easy but it is not, and cloud cover blows you away!

    Full Daylight Model Fitting

    Okay thats the new path, fit the full daylight model. The solver was a major coding hurdle but it is done now. If I feed the solver calculated LDR readings (quantitized to be fair) it has no problems solving latitude and longitude directly to within a few km (and the average cloud factor). A few hundred metres if the LDR reading are not quantized). If the data has consistent light reduction due to cloud cover, then no problem. It cloud cover is inconsistent (between morning and afternoon) then it will be a problem, you can't have everything! A better cloud cover model can be used but that is an optimisation for later.

    Although the TLS2591 would be ideal as it has higher resolution (16 bit) than the internal ADC (10 bit), the LDR and the internal ADC can be used. I calculate a 10k load resistor would be approprate (for my LDR). Here is the LDR ADC reading (red-right) versus Daylight (blue-left) Lux level used by my solver:


    The solver is based on the Nelder-Mead "ameoba method", but any steepest descent solver will work fine. Here is the inner core of the function (i.e. minmising Error2), the input variables to be optimised are CloudFactor, Lat and Lon:
      Error2=0.0;
      for (i=0;i<SampleSize;i++) {
        LT=SampleData[i][0]/60.0;
        SampleVoltage=(double)SampleData[i][2];
    
        // Determine GMST, RA and Dec
        GST=GMST(JD,LT-LTC);
        SunData(JD,LT-LTC,&RA,&Dec,&EoT);
    
        // Local Hour Angle and Altitude
        LHA=15.0*(GST-RA)+Lon;
        Alt=degrees(asin(sin(radians(Dec))*sin(radians(Lat))+cos(radians(Dec))*cos(radians(Lat))*cos(radians(LHA))));
        
        // Daylight intensity (Lux) and Sensor Voltage
        Lux=DayLight(Alt)*(1.0-CloudFactor);
        SensorVoltage=floor(1023.0/(1.0+10000.0/20000000.0*pow(Lux,0.75))+0.5);
    
        //Calculate Least Square Fit error
        Error2=Error2+(SensorVoltage-SampleVoltage)*(SensorVoltage-SampleVoltage);
      }

    Okay back to the Box code to include the new solver.

    Ported the code to the Box (Arduino) and a day of debugging. Using test data the Box reports a latitude within 30 km (longitude almost exact). No surprises that the latitude is off as it is the "hardest" to measure and calculate. Still using float, trials with double showed no improvement. The limitation is the "shallow bottom" of the search space and the effect of ADC quantisation.

    I am using the sum of absolute values rather than sum of squares because (1) the sum range is too great for float (at least initially so not really an issue) and (2) I want less weight to cloud and ADC quantisation noise.

    Interestingly enough if I use the square root of the absolute value then I get a very good result, within 1.3 km. I don't know if this is repeatable? Not convinced.

    Next is to rework the menu and sample storage code.

    Okay, the code has been rewritten and the Box is now collecting data. Hopefully no bugs.

    First Data

    The first data was off by 400 km. Not a great result. I am also getting some data corruption, perhaps it is time for a new RTC module. Anyway here is what the data looked like:

    The yellow line is the fitted model. In my mind it is a pretty good fit to the data (blue). The morning was cloudy and the afternoon clear. Also the sensor location is near a south-west facing window. The heel in the blue data would be direct sunlight on the sensor box in the afternoon. So okay, not the best box location. I could do much better.

    Two parameter SunFactor

    I reworked SunFactor for differences bewteen moring and afternoon sun light:

    The fit is quite good and the positional error is down to 120km. So really the next step is modifying the light sensor to be less sensitive to direct sunlight and putting the box in a better location.

    The solver is semi-stable. It needs...

    Read more »

  • Detecting cloud cover

    agp.cooper6 days ago 0 comments

    Detecting cloud cover

    As I have to use thresholds for dawn/dusk and/or sunrise/sunset then I need to know if the readings are affected by cloud cover.

    One option is to do record or integrate the maximum Lux for the day and then to adjust the sensor readings. Not ideal, and the LDR will need to be programmable over ranges (programmable pullup resistors).

    A quick search of the Internet and I found some people use of sky temperature (using far-infrared (heat) sensors. I have a couple of these sensor so it's a good option as it can be used at the dawk/dusk transition period. Calibrating the sensor to Lux attenuation is another project in itself.

    It also suggests that both the heat sensor and the Lux sensor should be pointed at the sky.

    I also came across the Adafruit TSL2591 Lux sensor. This seems to be a very good choice for replacing my primative LDR. One advantage of the Adafruit TSL2591 Lux sensor is that it has programmable Lux threshold interrupts. Now that sounds like a plan!

    Log graphs

    Well I learnt something about log graphs. They cannot be "trusted", well actually I cannot be trusted with them!

    Do you remember this one:

    The twilight transition really looks like a great place to do a measurement. Especially as the "standard" daylight Lux graph is almost identical:

    First lets put the "standard" in a linear graph:


    What happened to the twilight zone? Also it does not look right at 90 degrees, should it not flatten off? Well thats because it is a graph of Lux versus Altitude, not Time.
    Okay, lets look at a graph of Lux versus Time (for my place today) including sunrise, noon and sunset markers:
    Now that looks better.

    As we expect noisy data and fitting a least square fit model to the data is a way of dealing with the noise, what data points are important to measure? Answer, not twilight measurements! It would be the early morning and late afternoon measurements.
    The Lux/Time model is complex so we need to simplify. If we collect data for 60 minutes after sunrise or 60 minutes before sunset then simple quadratic has a very good fit:
    At sunrise or sunset the Lux is about 400 to 500 (456 based on my model).
    So after fitting a curve we can estimate sunrise (using a Lux of 450 and assuming clear conditions). And if conditions are not clear? Don't know that answer yet.
    If we can't be bothered with curve fitting then just use a threshold of 450.

    Magic

  • Useful references

    agp.cooper02/13/2018 at 05:44 0 comments

    Useful references

    I found this reference ("Variation of outdoor illumination as a function of solar elevation and light pollution") which has plots of dawn/dusk light intensity:

    https://www.nature.com/articles/srep26756#f2

    So there is a big difference between rural dawn and dusk (3 to 4 degrees!) but less for city dawn and dusk (about a degree). This will cause problems!

    Another reference (http://stjarnhimlen.se/comp/radfaq.html) had a sun Lux data set covering 90 degrees down to -18 degrees altitude! For -6 degrees (civil twilight) the light level is 3.4 Lux. This number would for a clear night. For the range -12 to 3 degrees I fitted the following cubic:

    Lux=10^(2.916+0.268942*Alt-0.03126*Alt^2-0.00156*Alt^3)

     So now I can map Time, Altitude, Lux and Sample Voltage:

    Notes:

    • the linear regression only uses sample voltages between 400 and 700
    • the assumed LDR dark resistance is 20M (I really should measure this!)
    • the LDR load resistor is 1.2M

    So two problems remain:

    • the dawn/dusk offset
    • non-ideal weather conditions

    No solution to the dawn/dusk offset except to use a constant (guess).

    Weather conditions

    A 50% reduction in light levels is shown below:

    This will cause problems as well.

    Oh Slope!

    I am coming to the conclusion that the slope is not that useful. In both cases the Lux level for the intercepts were the same (check the two graphs) so basically why bother! Again the slope like the inflection point is a function of the LDR maths. I should have picked this up when I was moving the curves left and right.

    Where to now?

    One option is to fit the entire daylight and twilight cycle. Entire curve will provide infromation on average/peak light levels. And I have a full model that I can fit.

    AlanX

  • Processing the Data

    agp.cooper02/11/2018 at 04:34 0 comments

    Processing the data

    Before rebuilding the Box for outside use (so the sensor has equal sensitivity in all directions), I throught I should try an tie down the data processing.

    What I want it to do:

    • Identify the transition curve pairs with a minimum of false signals.
    • Save the full transition curve pair windows for post analysis.
    • Consider noise reduction strategies.
    • Identify relationships between the transition curve and alititude.
    • Identify relationships that have the potential for "best" accuracy.
    • Save the results for later post analysis

    Transition curve identification and noise removal

    The image below shows the problem of artifical light on the signal an the potential for false signals:

    Properties of the transition curve:

    Twlight time:

    (source: http://individual.utoronto.ca/kalendis/twilight/MST.pdf)

    So the transition can take anywhere between 23 and 120 minutes for latitudes down to +/-60 degrees. So we could use a band pass filter between 360uHz and 69uHz (micro Hz!).

    Tried the bandpass, does not work very well. Used an exponential filter with a long period (a=1-1/60). With a 5 minute false trigger it works well:

    With a 70 minute window it works fine for my data. At a latitude of +/-60 degrees it will be far to short. Certainly possible to measure the curve slope at this time to adjust the adjust the window period for the next cycle.

    Okay, lets do that, here is a run where the initial window is 100 minutes with an adaptive code:

    The sensitively of the window to the false trigger is low and I am using a 5 minute window for the transit regression which seems fine.

    Here is the same with a initial window of 240 minutes:

    So I think I am halfway through my list of "wants".

    Defining Twilght

    Working on the second part of my wants list!

    Civil twilight or dawn/Dusk are defined as -6 degrees below the horizon. Great but it does not help my project. Sunrise/sunset could if I could set a sensor up to obsever the event (i.e. the direct light from a rising sun etc.) but that is not within the scope of my project.

    To date I have focused on the twilight "transition" but that needs refinement as it is about 30 minutes in duration. Around the transition time there is:

    • a night start/end point
    • a daylight start/end point
    • the mid-point
    • an inflection point, and
    • arbitary points in between.

    Forget about the "arbitary points in between" and the "mid-point" as they need to be calibrated for sensor setup AND the location.

    The points remaining that could be used are:

    • the night point
    • the day point
    • the inflection point

    So let us look at a typical sunrise:

    The sunrise area is a problem because the light level is high and the slope is quite flat. The curve is also not smooth? It will be difficult to get an accurate and consistent reading here. When you transpose the sensor sample voltage level into Lux, the inflection point does not exist:

    Even near dawn (the curve is a near perfect expontiental if you now where I am going):

    So the inflection point is due to the sensor response being non-linear. Here is the reverse sensor maths:

    • Lux=(20M*(1023-SampleReading)/1.2M/SampleReading)^(4/3)

    And to forward sensor maths:

    • Sample Voltage=1023/(1.2M/20M*(Lux)^(3/4) +1)

    In other words:

    • the LDR dark resistance is 20M (from the datasheet)
    • the LDR load resistance is 1.2M
    • the maximum sample reading is 1023
    • the LDR power law is 0.75 (from the data sheet).

    Now I have only two choices, daylight intercept or night-light intercept.

    The night-light intercept is interesting as it fits a logistic curve very well (the fitted red curve below):

    As an experiment I fitted an exponential Lux model AND the LDR model to the above data (red is the sunrise data and blue is the Lux/LDR model):

    Note: The tip of the pointers mark dawn and sunrise, and the mid-curve linear extrapolation...

    Read more »

  • You Got Me!

    agp.cooper02/10/2018 at 00:00 5 comments

    Arduino Variable Types

    I test a lot of my code in CodeBlocks using GCC on my Linux box.

    When done I copy and past into the Arduino sketch and just check that the out put matches before moving on.

    I am aware that arduino types are different:

    Arduino Linux

    Arduino Linux
    char    char
    byte    unsigned char
    int     short
    long    int
    

    Anyway the Arduino maths was not right!

    So checking various bits of my code I found the errors were too high.

    So as a check I printed a double calculation of effectively "4.45/24"  and to my horror I got "4.44999980". My "double" was a "float".

    Check the Arduino reference and guess what, "double" is "float"! Argh!

    Fortunately I can rearrange the key components of the calculations to reduce the code dependancy on double precision.

    Not Magic!

    Processing some data

    I collected two sunsets and a sun rise between 8 and 9 of February:

    The reversed altitudes for a threshold of 522 (precalculated for -3.5 degees from old data) was:

    • -4.48 degrees Dusk 8/2/2018
    • -5.67 degrees Dawn 9/2/2018
    • -4.55 degrees Dusk 9/2/2018

    So while the Dusk numbers appear repeatable, there is a difference between dawn and dusk light levels?

    Here is the dawn graph:

    Here is the dusk graph (aligned with dusk):

    Okay, I see, noise. I need to move down the curve (i.e. reduce the LDR load resistor). But the difference is only one time slot or 15 seconds. To be fair the Box uses a regression for the threshold rather than the raw data and in that case the difference would have been only 3 seconds. This is really a very good result but conditions were ideal.

    Dumb!

    Well I worked out why the sensor has a different reading for dawn versus dusk. It is mount vertically (on the inside) of a southern facing window and gets full morning sunlight. Adjacent to the sensor/window is an eastern facing wall that shields the sensor/window from all afternoon sunlight.

    Really, the sensor has to be mounted horizontally and above the roof so that light from all directions are treated equally.

    Anyway I can check the theory by placing the box on a western facing window.

    Now I probably have to separate the sensor from the box as the current box is not water proof and the idea of running cables for power etc. does not appeal to me.

    Magic

  • Test Runs

    agp.cooper02/05/2018 at 06:55 9 comments

      Test Runs

      I have a bug in my code. When I save my setup to the AT24C32 it overwrites an adjacent memory location (the sample count variable) that it should not. I need to find out why.

      Anyway I captured sunrise this morning and the data is quite useful to workout what to do next. It was a clear and moonless night. Practically perfect. Here is an early snapshot with twilight and sunrise marked:

      First it is pretty clear that I will not get good data if the sensor is not shielded from artifical light. It will swamp the area of most interest between civil twilight and sunrise.

      The area is of the most interest is because the curve is steep and relatively straight.

      I am running a 150k load resistor on my LDR sensor. Increasing this value will push down the curve increasing the relative range of the area of interest.

      Here is a second reading a few hours later:

      A Data Processing Plan

      1. Use a 90 minute window to find the transition. The average value of the window will be about 600 (based on the data above).
      2. Within that window calculate a linear least square fit with the data range trimmed to 950 (subject the background light levels including a full moon) and about 400 (based on the data above).

      You will get a line like this:

      The top intercept point (at 1023) looks like a pretty good point to use. It based on the steepest part of the curve and not on a threshold (i.e. in the graph above, the threshold for sunrise was about 600). A threshold is not a great idea because light levels change with weather conditions. The intercept was 10 minutes after twilight (or 20 minutes before sunrise). This implies to me I could get an accuracy in the order of +/-5 minutes (+/-10 km).

      AT24C32 Update

      I "patched" the code but the fault remains. Checked the library and the only thing I changed was to increase the delay to match the datasheet. I am busy collecting data so I cannot test the updated AT24C32 library. The library is not exactly compliant with the protocol (e.g. it does not listen for the data received handshake from the AT24C32) so the too short delay may be the issue. Works okay for single byte write but not "rapid fire".

      LDR Load Resistor

      After tonight's run I will decide if I want to increase the LDR load resistor. I think 1M will be about right.

      More Data Corruption

      I keep losing data! I lost my sunset. I will have to see if my modification to the library works. If not I will have to replace the RTC/AT2C32.

      Here is sunrise this morning. This time the sensor was placed looking out the window and shield from inside light:

      The noise at the bottom right of the graph is partial cloud cover.

      It really does look like that under ideal conditions the upper and lower intercept would pass through twilight and sunrise!

      The AT24C32 Fixed

      Increasing the delay from 5 ms to 10 ms (as per the datasheet). It did not work. I moved the settings and the sample counter to different EEPROM pages. I will likely have to but a new RTC.

      A Sunset run

      Increased the LDR load resistor to 1.2M (it was all I had). For this run I left the LCD LED on by accident, however back calculating the LDR resistance was 10M. The minimum dark resistance for the LDE is 20M (minimum), so the dark level might increase from 915 to 965 under ideal conditions:

      Still the fitted line is within 3 minutes of civil sunset using the 965 intercept.

      I quite like this set up as the daylight noise is practically gone. Replacing the 1.2M with a 1M resistor will increase the dark level...

    Read more »

  • The Box

    agp.cooper02/04/2018 at 11:37 0 comments

    The Box

    Need a box to put the electronics in. I have a DeltaCAD macro that generates boxes for a laser cutter. For this project I modified the macro for a router type CNC machine. Makes sense as I don't have a laser cutter but I do have a router CNC machine.

    The Macro (CaseII.bas)

    The macro is not very clever. It is just a long list of line commands with each point carefully calculated and checked. Here is the macro:

    Option Explicit
    
    ' Maximize the window, close any existing drawing without saving and start a new drawing.
    dcSetDrawingWindowMode dcMaximizeWin
    dcCloseWithoutSaving
    dcNew ""
    
    ' Drawing units and scale
    dcSetDrawingScale 1
    dcSetDrawingUnits dcMillimeters
    
    Call Main
    End
    
    
    Sub Main
    
    ' Units are mm
    Dim BoxSheetThickness As Double
    Dim BoxFlangeWidth As Double
    Dim FlangeHoleDiameter As Double
    Dim OverCutWidth As Double
    Dim EndMillDiameter As Double
    
    ' Inside dimensions
    Dim BoxDepth As Double
    Dim BoxWidth As Double
    Dim BoxHeight As Double
    
    BoxDepth=25.0
    BoxWidth=25.0
    BoxHeight=24.0
    FlangeHoleDiameter=3.0
    BoxSheetThickness=3.0
    BoxFlangeWidth=2*BoxSheetThickness+FlangeHoleDiameter
    OverCutWidth=0.0
    EndMillDiameter=BoxSheetThickness/2
    
    Begin Dialog DesignInput 50,50,160,110, "Input data for the Box"
    
      Text 10,10,80,10, "Inside box depth"
      Text 10,20,80,10, "Inside box width"
      Text 10,30,80,10, "Inside box height"
      Text 10,40,80,10, "Box sheet thickness"
      Text 10,50,80,10, "Box flange width"
      Text 10,60,80,10, "Flange hole diameter"
      Text 10,70,80,10, "Over (under) cut width"
      Text 10,80,80,10, "End Mill Diameter"
      TextBox 90,10,45,10, .bd
      TextBox 90,20,45,10, .bw
      TextBox 90,30,45,10, .bh
      TextBox 90,40,45,10, .st
      TextBox 90,50,45,10, .fw
      TextBox 90,60,45,10, .hd
      TextBox 90,70,45,10, .cw
      TextBox 90,80,45,10, .ed
      Text 140,10,15,10, "mm"
      Text 140,20,15,10, "mm"
      Text 140,30,15,10, "mm"
      Text 140,40,15,10, "mm"
      Text 140,50,15,10, "mm"
      Text 140,60,15,10, "mm"
      Text 140,70,15,10, "mm"
      Text 140,80,15,10, "mm"
    
      OKButton 10,95,30,10
      CancelButton 90,95,30,10
    
    End Dialog
    
    'Initialize 
    Dim DesignData As DesignInput
    Dim button As Long
    
    DesignData.bd=BoxDepth
    DesignData.bw=BoxWidth
    DesignData.bh=BoxHeight
    DesignData.st=BoxSheetThickness
    DesignData.fw=BoxFlangeWidth
    DesignData.hd=FlangeHoleDiameter
    DesignData.cw=OverCutWidth
    DesignData.ed=EndMillDiameter
    
    button=Dialog(DesignData)
    
    If (button=True) Then 
    
      BoxDepth=CDbl(DesignData.bd)
      BoxWidth=CDbl(DesignData.bw)
      BoxHeight=CDbl(DesignData.bh)
      BoxSheetThickness = CDbl(DesignData.st)
      BoxFlangeWidth=CDbl(DesignData.fw)
      FlangeHoleDiameter=CDbl(DesignData.hd)
      OverCutWidth=CDbl(DesignData.cw)
      EndMillDiameter=CDbl(DesignData.ed)
      
      Dim i As Long
      Dim XOffset As Double
      Dim YOffset As Double
    
      dcSetLineParms dcBLACK, dcSOLID, dcTHIN
      dcSetCircleParms dcBLACK, dcSOLID, dcTHIN
    
      XOffset=5.0  
      YOffset=5.0  
    
      ' Top and bottom
      For i=1 to 2
    
        dcSetCircleParms dcBLUE, dcSOLID, dcTHIN
    
        ' Outside (outline)
        dcCreateLine _
          XOffset-OverCutWidth/2, _
          YOffset-OverCutWidth/2+BoxSheetThickness/2, _
          XOffset-OverCutWidth/2+BoxSheetThickness/2, _
          YOffset-OverCutWidth/2
    
        dcCreateLine _
          XOffset-OverCutWidth/2+BoxSheetThickness/2, _
          YOffset-OverCutWidth/2, _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth-BoxSheetThickness/2, _
          YOffset-OverCutWidth/2
    
        dcCreateLine _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth-BoxSheetThickness/2, _
          YOffset-OverCutWidth/2, _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth, _
          YOffset-OverCutWidth/2+BoxSheetThickness/2
    
        dcCreateLine _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth, _
          YOffset-OverCutWidth/2+BoxSheetThickness/2, _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth, _
          YOffset+OverCutWidth/2+BoxDepth+2*BoxFlangeWidth-BoxSheetThickness/2
    
        dcCreateLine _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth, _
          YOffset+OverCutWidth/2+BoxDepth+2*BoxFlangeWidth-BoxSheetThickness/2, _
          XOffset+OverCutWidth/2+BoxWidth+2*BoxFlangeWidth-BoxSheetThickness/2, _
          YOffset+OverCutWidth/2+BoxDepth+2*BoxFlangeWidth
    
        dcCreateLine _
     XOffset+OverCutWidth/...
    Read more »

  • The Mathematics

    agp.cooper02/04/2018 at 09:04 0 comments

      Julian Day (JD) number

      Most astronomy based formulas use the Julian day number as an input.

      I use the a modified 2000 Julian Day number, the number of days since mid-night of the 1st of January 2000:

      double J2000(long Year,long Month,long Day,double UT) {
        double JD;
        
        // Julian (2000) Day Number
        JD=(367*Year-(7*(Year+(Month+9)/12))/4+275*Month/9+Day);
        JD=JD-730531.5;
        JD=JD+UT/24.0;
        
        return JD;
      }
      

      The formular is correct even for dated before 1/1/2000.

      One of the inputs to the function is UT or Universal Time. UT is the time in Greenwich UK. I will write more about this later.

      Note that Year, Month and Day are type long. Therefore integer division rules apply.

      Reverse of the Julian Day number

      I found a source for calculating the reverse of the 2000 Julian Day Number:

      Source: http://mathforum.org/library/drmath/view/51907.htm

        unsigned long JD,p,q,r,s,t,u,v;
      
        p=JD+2451545+68569;
        q=4*p/146097;
        r=p-(146097*q+3)/4;
        s=4000*(r+1)/1461001;
        t=r-1461*s/4+31;
        u=80*t/2447;
        v=u/11;
        Year=100*(q-49)+s+v;
        Month=u+2-12*v;
        Day=t-2447*u/80;

      It may be useful when checking/calculating valid dates. I checked the code and it good from 1/1/2000 to 31/12/2099.

      Greenwich Mean Sideral Time (GMST)

      This is the time of the stars! You would be aware that there are 365.25 days in a solar year. But the stars rotate around the earth 366.25 times a year. The extra rotation is because the earth also rotates around the sun as well as rotating around it axis. Here is the GMST function:

      double GMST(double JD,double UT) {
        double gmst;
        
        // Greenwich Mean Sidereal Time (~36625/36525)
        gmst=18.697374558+24.06570982441908*(JD+UT/24.0);
        gmst=fmod(fmod(gmst,24.0)+24.0,24.0);
        
        return gmst;
      }

       We need GMST even for the sun as the sun is cosidered just another start. But the Right Ascension and Declination moves a lot more than ordinary stars.

      Sun Right Ascension (RA) and Declination (Dec)

      Now we can calculate the position of the sun amongst the stars. I use a series of formulas from my copy of the "The Astronomical Amanac (1986)" but I have seen the same formulas on the Internet:

      void SunData(double JD,double *RA,double *Dec,double *EoT) {
        double L,g,a,b,c;
        
        /* Source: Page C24 from "The Astronomical Almanac" (1986) */
        // Mean longitude of Sun, corrected for aberration
        L=280.460+0.9856474*JD;
        L=fmod(fmod(L,360.0)+360.0,360.0);
        // Mean anomaly
        g=357.528+0.9856003*JD;
        g=fmod(fmod(g,360.0)+360.0,360.0);
        // Ecliptic longitude
        a=L+1.915*sin(g*M_PI/180)+0.020*sin(2*g*M_PI/180);
        a=fmod(fmod(a,360.0)+360.0,360.0);
        // Ecliptic latitude
        b=0;
        // Obliquity of ecliptic
        c=23.439-0.0000004*JD;
        c=fmod(fmod(c,360.0)+360.0,360.0);
        // Right ascension
        (*RA)=atan2(cos(c*M_PI/180)*sin(a*M_PI/180),cos(a*M_PI/180))*12/M_PI;
        if ((*RA)<0) (*RA)=(*RA)+24.0;
        // Equation of Time
        (*EoT)=(L/15.0-(*RA));
        // Declination
        (*Dec)=asin(sin(c*M_PI/180)*sin(a*M_PI/180))*180/M_PI;
      }

      The procedure accepts modified Julian Day number (JD) and returns Right Ascension (RA), Declination (Dec) and the Equarion of Time (EoT). The units for RA and EoT are hours but the units for Declination is degrees.

      Checking the formulas

      It is really easy to make a minor mistake with these formulas so it is important to check the results against a worked example. The Internet has many calculators that you can check your results against.

      Sunrise and sunset

      Sunrise and sunset times are published on the Internet for any place that you live.

      On the 20th of January 2018 in Perth Western Australia the Sunrise was 5:30 and Sunset was 19:24. Local noon was therefore 12:27. The local time is 8 hours ahead of UT.

      Calculating  longitude

      The steps are:

      1.   SunRise=5.50-8;
      2.  ...
    Read more »

View all 8 project logs

Enjoy this project?

Share

Discussions

ris wrote 02/12/2018 at 19:36 point

This project has got me thinking a lot. I wonder how many different global radio clock sources (https://en.wikipedia.org/wiki/Radio_clock#List_of_radio_time_signal_stations) someone could get a low power microcontroller to recognize, with minimal, cheap external components. And thus how much global coverage could be achieved.

  Are you sure? yes | no

agp.cooper wrote 8 hours ago point

Sorry, I missed your comment. The radio clock is an interresting option and I have looked at it but for me (in Australia) it is not an option (no radio-time transmitters). If your are in EU, US, Japan or China then it will work. At 60kHz the range is about 2,000km. One guy "cheated' by using his PC to get the time (from the internet) and transmitting his own 60kHz signal! Local AM radio stations used to transmit pips on the hour but not anymore.
A RTC is really the only option.
AlanX

  Are you sure? yes | no

Ted Yapo wrote 7 hours ago point

You can get accurate time from GPS ... oh, wait :-)

  Are you sure? yes | no

Mike Szczys wrote 02/05/2018 at 17:01 point

Your tests of artificial light interference are really interesting! I also appreciated you walking us through the equations you've implemented.

@jaromir.sukuba it's great to see others trying out the geolocator techniques you were working on!

  Are you sure? yes | no

Similar Projects

Does this project spark your interest?

Become a member to follow this project and never miss any updates