Jump to content

standalone moon phase calculation


Gianni
 Share

Recommended Posts

I'm trying to translate this C listing found here http://www.voidware.com/phase.c to AutoIt.
.... i'm working on debugging a messy translation done by me, but it still doesn't work...
My question is: how this snippet could be translated to AutoIt??

double moon_phase(int year,int month,int day, double hour, int* ip)
{
    /*
      Calculates more accurately than Moon_phase , the phase of the moon at
      the given epoch.
      returns the moon phase as a real number (0-1)
      */

    double j= Julian(year,month,(double)day+hour/24.0)-2444238.5;
    double ls = sun_position(j);
    double lm = moon_position(j, ls);

    double t = lm - ls;
    if (t < 0) t += 360;
    *ip = (int)((t + 22.5)/45) & 0x7;
    return (1.0 - cos((lm - ls)*RAD))/2;
}

in particular this line :

*ip = (int)((t + 22.5)/45) & 0x7;

Thanks

Edited by Chimp

 

image.jpeg.9f1a974c98e9f77d824b358729b089b0.jpeg Chimp

small minds discuss people average minds discuss events great minds discuss ideas.... and use AutoIt....

Link to comment
Share on other sites

I'm a C noop and I would interpret it this way:

Func Moon_Phase($year, $month, $day, $hour, ByRef $ip)
    Local $j = Julian($year, $month, $day + $hour / 24.0) - 2444238.5
    Local $ls = Sun_Position($j)
    Local $lm = Moon_Position($j, $ls)
    Local $t = $lm - $ls
    If ($t < 0) Then $t += 360
    $ip = BitAND(Int(($t + 22.5) / 45), 0x7)
    Return (1.0 - Cos(($lm - $ls) * $fRad)) / 2
EndFunc   ;==>Moon_Phase

I also tried to translate it to AutoIt but I don't know what "->" means:

now->day = c - e - g1;
    now->hour = (jd - jdi) * 24.0;
    if (g <= 13) now->month = g - 1;
    else now->month = g - 13;
    if (now->month > 2) now->year = d - 4716;
    else now->year = d - 4715;

 

Please don't send me any personal message and ask for support! I will not reply!

Selection of finest graphical examples at Codepen.io

The own fart smells best!
Her 'sikim hıyar' diyene bir avuç tuz alıp koşma!
¯\_(ツ)_/¯  ٩(●̮̮̃•̃)۶ ٩(-̮̮̃-̃)۶ૐ

Link to comment
Share on other sites

Hi @UEZ, thanks for your translation, it seems that it's ok that way.
Regarding the "->" syntax, seems that is a way to access variables within a structure.
I tried using a global array instead, and seems that it works as well.

The translation now is complete...(?), but according to the result produced by the original executable found herehttp://www.voidware.com/moon_phase.htm (see phase.exe near the bottom of the page) the results by this listing, are slight different.
I don't know if that is due to some imprecise in my translation or different precision in AutoIt math.... (???)

Anyway the moon phases provided by this listing are correct, have a look here for example to check: http://cycletourist.com/moon/lunar_calendar.php

fixes or improvements are welcome

Thanks everybody

EDIT: Listing corrected as suggested by LarsJ in post #4 and now works OK. (Thanks LarsJ :))

; http://www.voidware.com/phase.c
; /* standalone moon phase calculation */

Const $PI = 3.1415926535897932384626433832795
Const $RAD = ($PI / 180.0)
Const $SMALL_FLOAT = (1E-12)

Global $TimePlace[4] ; [0] Year ; [1] Month ; [2] ; Day ; [3] Hour
Global $ip

main()

Func JulianToDate($jd)

    $jd += 0.5
    $jdi = Int($jd)
    If ($jdi > 2299160) Then
        $a = Int(($jdi - 1867216.25) / 36524.25)
        $b = Int($jdi + 1 + $a - ($a / 4))
    Else
        $b = $jdi
    EndIf
    $c = $b + 1524;
    $d = Int(($c - 122.1) / 365.25) ;
    $e = Int(365.25 * $d) ;
    $g = Int(($c - $e) / 30.6001);
    $g1 = Int(30.6001 * $g) ;
    $TimePlace[2] = Int($c - $e - $g1) ; Day
    $TimePlace[3] = ($jd - $jdi) * 24 ; Hour
    If ($g <= 13) Then
        $TimePlace[1] = Int($g - 1) ; Month
    Else
        $TimePlace[1] = Int($g - 13) ; Month
    EndIf
    If $TimePlace[1] > 2 Then ; Month
        $TimePlace[0] = Int($d - 4716)
    Else
        $TimePlace[0] = Int($d - 4715); Year
    EndIf
EndFunc   ;==>JulianToDate

; double
Func Julian($year, $month, $day)

    Local $y = $year, $m = $month, $d = $day

    #cs   /*
        Returns the number of julian days for the specified day.
    #ce   */

    If ($month < 3) Then
        $year -= 1
        $month += 12
    EndIf

    If $year > 1582 Or ($year = 1582 And $month > 10) Or ($year = 1582 And $month = 10 And $day > 15) Then
        $a = Int($year / 100) ;
        $b = 2 - $a + Int($a / 4)
    EndIf

    $c = Int(365.25 * $year);
    $e = Int(30.6001 * ($month + 1))
    Return $b + $c + $e + $day + 1720994.5

EndFunc   ;==>Julian

Func sun_position($j)

    Local $dl
    $n = 360 / 365.2422 * $j
    $i = Int($n / 360)
    $n = $n - $i * 360.0
    $x = $n - 3.762863 ;
    If ($x < 0) Then $x += 360
    $x *= $RAD
    $e = $x

    #cs ; This while - wend replaced by the following Do Until loop
        ; as suggested by @LarsJ in post #4
    While Abs($dl) >= $SMALL_FLOAT
        $dl = $e - 0.016718 * Sin($e) - $x
        $e = $e - $dl / (1 - 0.016718 * Cos($e))
    WEnd
    #ce
    
    Do
        $dl = $e - 0.016718 * Sin($e) - $x
        $e = $e - $dl / (1 - 0.016718 * Cos($e))
    Until Abs($dl) < $SMALL_FLOAT

    $v = 360 / $PI * ATan(1.01686011182 * Tan($e / 2))
    $l = $v + 282.596403
    $i = Int($l / 360)
    $l = $l - $i * 360.0
    Return $l
EndFunc   ;==>sun_position

Func moon_position($j, $ls)
    ; /* ls = sun_position(j) */
    $ms = 0.985647332099 * $j - 3.762863
    If ($ms < 0) Then $ms += 360.0
    $l = 13.176396 * $j + 64.975464
    $i = Int($l / 360)
    $l = $l - $i * 360.0
    If ($l < 0) Then $l += 360.0
    $mm = $l - 0.1114041 * $j - 349.383063
    $i = Int($mm / 360)
    $mm -= $i * 360.0
    $n = 151.950429 - 0.0529539 * $j
    $i = Int($n / 360)
    $n -= $i * 360.0
    $ev = 1.2739 * Sin((2 * ($l - $ls) - $mm) * $RAD)
    $sms = Sin($ms * $RAD)
    $ae = 0.1858 * $sms
    $mm += $ev - $ae - 0.37 * $sms
    $ec = 6.2886 * Sin($mm * $RAD)
    $l += $ev + $ec - $ae + 0.214 * Sin(2 * $mm * $RAD)
    $l = 0.6583 * Sin(2 * ($l - $ls) * $RAD) + $l
    Return $l
EndFunc   ;==>moon_position

Func moon_phase($year, $month, $day, $hour, ByRef $ip)
    #cs  /*
        Calculates more accurately than Moon_phase , the phase of the moon at
        the given epoch.
        returns the moon phase as a real number (0-1)
    #ce  */

    $j = Julian($year, $month, $day + $hour / 24.0) - 2444238.5
    $ls = sun_position($j)
    $lm = moon_position($j, $ls)

    $t = $lm - $ls
    If ($t < 0) Then $t += 360 ;
    ; * ip = (int) ((t + 22.5) / 45) & 0x7 ;  <----- ????
    $ip = BitAND(Int(($t + 22.5) / 45), 0x7)
    Return (1.0 - Cos(($lm - $ls) * $RAD)) / 2
EndFunc   ;==>moon_phase

Func nextDay($y, $m, $d, $dd)
    $jd = Julian($y, $m, $d)
    $jd += $dd
    JulianToDate($jd) ; set result in $TimePlace[]
    #cs
        $y = $TimePlace[0] ; Year
        $m = $TimePlace[1] ; month
        $d = $TimePlace[2] ; day
    #ce
EndFunc   ;==>nextDay

Func main() ; {
    $y = @YEAR
    $m = @MON
    $d = @MDAY

    $step = 1
    $begun = 0

    $pmax = 0
    $pmin = 1
    $plast = 0

    ConsoleWrite("tabulation of the phase of the moon for one month" & @CRLF & @CRLF)

    ConsoleWrite("Year: " & $y & @CRLF)
    ConsoleWrite("month: " & $m & @CRLF)

    $d = 1
    $m0 = Int($m)

    ConsoleWrite(@CRLF & "Date       Time   Phase Segment" & @CRLF)

    While True
        For $h = 0 To 23 Step $step

            $p = moon_phase($y, $m, $d, $h, $ip)

            If ($begun) Then
                If ($p > $plast And $p > $pmax) Then
                    $pmax = $p
                    $ymax = $y
                    $mmax = $m
                    $dmax = $d
                    $hmax = $h
                ElseIf ($pmax) Then ; {
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00          (fullest)\n", $ymax, $mmax, $dmax, $hmax))
                    $pmax = 0
                EndIf

                If ($p < $plast And $p < $pmin) Then
                    $pmin = $p
                    $ymin = $y
                    $mmin = $m
                    $dmin = $d
                    $hmin = $h
                ElseIf ($pmin < 1) Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00          (newest)\n", $ymin, $mmin, $dmin, $hmin));
                    $pmin = 1.0
                EndIf

                If ($h = 16) Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00 %5.1f%%   (%d)\n", $y, $m, $d, $h, Floor($p * 1000 + 0.5) / 10, $ip))
                EndIf
            Else
                $begun = 1
            EndIf
            $plast = $p;
        Next
        nextDay($y, $m, $d, 1.0);

        $y = $TimePlace[0] ; Year
        $m = $TimePlace[1] ; month
        $d = $TimePlace[2] ; day
        If ($m <> $m0) Then ExitLoop ; break
    WEnd
    Return 0;
EndFunc   ;==>main

 

Edited by Chimp

 

image.jpeg.9f1a974c98e9f77d824b358729b089b0.jpeg Chimp

small minds discuss people average minds discuss events great minds discuss ideas.... and use AutoIt....

Link to comment
Share on other sites

Replace the While-WEnd loop in sun_position()

While Abs($dl) >= $SMALL_FLOAT
    $dl = $e - .016718 * Sin($e) - $x
    $e = $e - $dl / (1 - .016718 * Cos($e))
WEnd

with a Do-Until loop:

Do
    $dl = $e - .016718 * Sin($e) - $x
    $e = $e - $dl / (1 - .016718 * Cos($e))
Until Abs($dl) >= $SMALL_FLOAT

It guarantees that the calculation in the loop is done at least once.

 

On second thought: Just drop the loop and make sure that the calculation is done once.

Edited by LarsJ
Link to comment
Share on other sites

Shouldn't it be

Do
    $dl = $e - .016718 * Sin($e) - $x
    $e = $e - $dl / (1 - .016718 * Cos($e))
Until Abs($dl) < $SMALL_FLOAT

instead?

 

Btw, this is what I did yesterday but it's not working:

Global Const $fPI = ACos(-1), $fRad = $fPI / 180, $fSmall_Float = 1E-12
Global $tTimePlace = DllStructCreate("struct;int year;int month;int day;double hour;endstruct")

Func JulianToDate(ByRef $tTimePlace, $jd)
    Local $jdi, $a, $b, $c, $d, $e, $g, $g1
    $jd += 0.5
    $jdi = $jd
    If ($jdi > 2299160.0) Then
        $a = ($jdi - 1867216.25) / 36524.25
        $b = $jdi + 1 + $a - $a / 4
    Else
        $b = $jdi
    EndIf
    $c = $b + 1524
    $d = ($c - 122.1) / 365.25
    $e = 365 * $d
    $g = ($c - $e) / 30.6001
    $g1 = 30.6001 * $g

    $tTimePlace.day = Int($c - $e - $g1)
    $tTimePlace.hour = ($jd - $jdi) * 24.0
    If $g <= 13 Then
        $tTimePlace.month = Int($g - 1)
    Else
        $tTimePlace.month = Int($g - 13)
    EndIf
    If $tTimePlace.month > 2 Then
        $tTimePlace.year = Int($d - 4716)
    Else
        $tTimePlace.year = Int($d - 4715)
    EndIf
EndFunc   ;==>JulianToDate

Func Julian($year, $month, $day)
    Local $a, $b, $c, $e
    If $month < 3 Then
        $year -= 1
        $month += 12
    EndIf
    If ($year > 1582 Or ($year = 1582 And $month > 10) Or ($year = 1582 And $month = 10 And $day > 15)) Then
        $a = Int($year / 100)
        $b = Int(2 - $a + $a / 4)
    EndIf
    $c = Int(365.25 * $year)
    $e = Int(30.6001 * ($month + 1))
    Return $b + $c + $e + $day + 1720994.5
EndFunc   ;==>Julian

Func Sun_Position($j)
    Local $n, $x, $e, $l, $dl, $v, $m2, $i
    $n = 360 / 365.2422 * $j
    $i = Int($n / 360)
    $n = $n - $i * 360
    $x = $n - 3.762863
    If $x < 0 Then $x += 360
    $x *= $fRad
    $e = $x
    Do
        $dl = $e - 0.016718 * Sin($e) - $x
        $e = $e - $dl / (1 - 0.016718 * Cos($e))
    Until Abs($dl) < $fSmall_Float
    $v = 360 / $fPI * ATan(1.01686011182 * Tan($e / 2))
    $l = $v + 282.596403
    $i = Int($l / 360)
    $l = $l - $i * 360.0
    Return $l
EndFunc   ;==>Sun_Position

Func Moon_Position($j, $ls)
    Local $ms, $l, $mm, $n, $ev, $sms, $z, $x, $lm, $bm, $ae, $ec, $d, $ds, $as, $dm, $i
    $ms = 0.985647332099 * $j - 3.762863
    If ($ms < 0) Then $ms += 360.0
    $l = 13.176396 * $j + 64.975464
    $i = Int($l / 360)
    $l = $l - $i * 360.0
    If ($l < 0) Then $l += 360.0
    $mm = $l - 0.1114041 * $j - 349.383063
    $i = Int($mm / 360)
    $mm -= $i * 360.0
    $n = 151.950429 - 0.0529539 * $j
    $i = Int($n / 360)
    $n -= $i * 360.0
    $ev = 1.2739 * Sin((2 * ($l - $ls) - $mm) * $fRad)
    $sms = Sin($ms * $fRad)
    $ae = 0.1858 * $sms
    $mm += $ev - $ae - 0.37 * $sms
    $ec = 6.2886 * Sin($mm * $fRad)
    $l += $ev + $ec - $ae + 0.214 * Sin(2 * $mm * $fRad)
    $l = 0.6583 * Sin(2 * ($l - $ls) * $fRad) + $l
    Return $l
EndFunc   ;==>Moon_Position

Func Moon_Phase($year, $month, $day, $hour, ByRef $ip)
    Local $j = Julian($year, $month, $day + $hour / 24.0) - 2444238.5
    Local $ls = Sun_Position($j)
    Local $lm = Moon_Position($j, $ls)
    Local $t = $lm - $ls
    If ($t < 0) Then $t += 360
    $ip = BitAND(Int(($t + 22.5) / 45), 0x7)
    Return (1.0 - Cos(($lm - $ls) * $fRad)) / 2
EndFunc   ;==>Moon_Phase

Func NextDay(ByRef $y, ByRef $m, ByRef $d, $dd)
    Local $jd = Julian($y, $m, $d)
    $jd += $dd
    JulianToDate($tTimePlace, $jd)
    $y = $tTimePlace.year
    $m = $tTimePlace.month
    $d = $tTimePlace.day
EndFunc   ;==>NextDay

Func Main($y = @YEAR, $m = @MON)
    Local $d, $m0, $h, $i, $p, $ip = 0
    Local $step = 1, $begun = 0, $pmax = 0, $pmin = 1, $plast = 0
    Local $ymax, $mmax, $dmax, $hmax, $ymin, $mmin, $dmin, $hmin
    $d = 1
    $m0 = $m
    ConsoleWrite("Date       Time   Phase Segment" & @CRLF)
    Do
        $ip = 0
        For $h = 0 To 23 Step $step
            $p = Moon_Phase($y, $m, $d, $h, $ip)
            If $begun Then
                If ($p > $plast And $p > $pmax) Then
                    $pmax = $p
                    $ymax = Int($y)
                    $mmax = Int($m)
                    $dmax = Int($d)
                    $hmax = Int($h)
                ElseIf ($pmax) Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00          (fullest)\n", $ymax, $mmax, $dmax, $hmax))
                    $pmax = 0
                ElseIf $p < $plast And $p < $pmin Then
                    $pmin = $p
                    $ymin = Int($y)
                    $mmin = Int($m)
                    $dmin = Int($d)
                    $hmin = Int($h)
                ElseIf $pmin < 1 Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00          (newest)\n", $ymin, $mmin, $dmin, $hmin))
                    $pmin = 1.0
                EndIf
                If $h = 16 Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00 %5.1f%%   (%d)\n", $y, $m, $d, $h, Floor($p * 1000 + 0.5) / 10, $ip))
                EndIf
            Else
                $begun = 1
            EndIf
            $plast = $p
        Next
        NextDay($y, $m, $d, 1.0)
        If ($m <> $m0) Then ExitLoop
    Until False
EndFunc   ;==>Main


Main()

I cannot find yet the mistakes...

Edited by UEZ

Please don't send me any personal message and ask for support! I will not reply!

Selection of finest graphical examples at Codepen.io

The own fart smells best!
Her 'sikim hıyar' diyene bir avuç tuz alıp koşma!
¯\_(ツ)_/¯  ٩(●̮̮̃•̃)۶ ٩(-̮̮̃-̃)۶ૐ

Link to comment
Share on other sites

UEZ, You are right about the loop.

 

Your calculations in JulianToDate

$jd += 0.5
$jdi = $jd
If ($jdi > 2299160.0) Then
    $a = ($jdi - 1867216.25) / 36524.25
    $b = $jdi + 1 + $a - $a / 4
Else
    $b = $jdi
EndIf
$c = $b + 1524
$d = ($c - 122.1) / 365.25
$e = 365 * $d
$g = ($c - $e) / 30.6001
$g1 = 30.6001 * $g

will definitely not give the same results as these calculations by Chimp

$jd += 0.5
$jdi = Int($jd)
If ($jdi > 2299160) Then
    $a = Int(($jdi - 1867216.25) / 36524.25)
    $b = Int($jdi + 1 + $a - ($a / 4))
Else
    $b = $jdi
EndIf
$c = $b + 1524;
$d = Int(($c - 122.1) / 365.25) ;
$e = Int(365.25 * $d) ;
$g = Int(($c - $e) / 30.6001);
$g1 = Int(30.6001 * $g) ;

 

Link to comment
Share on other sites

Indeed LarsJ. :thumbsup:

When I run my script with the modification from above it seems to be "more or less equal" as from the C code:

AutoIt:

Date       Time   Phase Segment
2016/05/01 16:00  33.7%   (6)
2016/05/02 16:00  23.4%   (7)
2016/05/03 16:00  14.2%   (7)
2016/05/04 16:00   6.9%   (7)
2016/05/05 16:00   2.0%   (0)
2016/05/06 16:00   0.0%   (0)
2016/05/07 16:00   1.1%   (0)
2016/05/08 16:00   5.1%   (1)
2016/05/09 16:00  11.6%   (1)
2016/05/10 16:00  19.8%   (1)
2016/05/11 16:00  29.2%   (1)
2016/05/12 16:00  39.3%   (2)
2016/05/13 16:00  49.4%   (2)
2016/05/14 16:00  59.3%   (2)
2016/05/15 16:00  68.7%   (2)
2016/05/16 16:00  77.2%   (3)
2016/05/17 16:00  84.6%   (3)
2016/05/18 16:00  90.8%   (3)
2016/05/19 16:00  95.5%   (3)
2016/05/20 16:00  98.6%   (4)
2016/05/21 16:00 100.0%   (4)
2016/05/21 21:00          (fullest)
2016/05/06 19:00          (newest)
2016/05/22 16:00  99.4%   (4)
2016/05/23 16:00  97.0%   (4)
2016/05/24 16:00  92.6%   (5)
2016/05/25 16:00  86.5%   (5)
2016/05/26 16:00  78.7%   (5)
2016/05/27 16:00  69.5%   (5)
2016/05/28 16:00  59.2%   (6)
2016/05/29 16:00  48.3%   (6)
2016/05/30 16:00  37.1%   (6)
2016/05/31 16:00  26.4%   (7)

C:

tabulation of the phase of the moon for one month

year: 2016
month: 5

Date       Time   Phase Segment
2016/05/01 16:00  33.7%   (6)
2016/05/02 16:00  23.4%   (7)
2016/05/03 16:00  14.2%   (7)
2016/05/04 16:00   6.9%   (7)
2016/05/05 16:00   2.0%   (0)
2016/05/06 16:00   0.0%   (0)
2016/05/06 19:00          (newest)
2016/05/07 16:00   1.1%   (0)
2016/05/08 16:00   5.1%   (1)
2016/05/09 16:00  11.6%   (1)
2016/05/10 16:00  19.8%   (1)
2016/05/11 16:00  29.2%   (1)
2016/05/12 16:00  39.3%   (2)
2016/05/13 16:00  49.4%   (2)
2016/05/14 16:00  59.3%   (2)
2016/05/15 16:00  68.7%   (2)
2016/05/16 16:00  77.2%   (3)
2016/05/17 16:00  84.6%   (3)
2016/05/18 16:00  90.8%   (3)
2016/05/19 16:00  95.5%   (3)
2016/05/20 16:00  98.6%   (4)
2016/05/21 16:00 100.0%   (4)
2016/05/21 21:00          (fullest)
2016/05/22 16:00  99.4%   (4)
2016/05/23 16:00  97.0%   (4)
2016/05/24 16:00  92.6%   (5)
2016/05/25 16:00  86.5%   (5)
2016/05/26 16:00  78.7%   (5)
2016/05/27 16:00  69.5%   (5)
2016/05/28 16:00  59.2%   (6)
2016/05/29 16:00  48.3%   (6)
2016/05/30 16:00  37.1%   (6)
2016/05/31 16:00  26.4%   (7)

Time and sequment is not fully ok...

Edited by UEZ

Please don't send me any personal message and ask for support! I will not reply!

Selection of finest graphical examples at Codepen.io

The own fart smells best!
Her 'sikim hıyar' diyene bir avuç tuz alıp koşma!
¯\_(ツ)_/¯  ٩(●̮̮̃•̃)۶ ٩(-̮̮̃-̃)۶ૐ

Link to comment
Share on other sites

Link to comment
Share on other sites

Ok, this seems to be creating the same result as the C code:

Global Const $fPI = ACos(-1), $fRad = $fPI / 180, $fSmall_Float = 1E-12
Global $tTimePlace = DllStructCreate("struct;int year;int month;int day;double hour;endstruct")

Func JulianToDate(ByRef $tTimePlace, $jd)
    Local $jdi, $a, $b, $c, $d, $e, $g, $g1
    $jd += 0.5
    $jdi = Int($jd)
    If ($jdi > 2299160) Then
        $a = Int(($jdi - 1867216.25) / 36524.25)
        $b = Int($jdi + 1 + $a - ($a / 4))
    Else
        $b = $jdi
    EndIf
    $c = $b + 1524;
    $d = Int(($c - 122.1) / 365.25)
    $e = Int(365.25 * $d)
    $g = Int(($c - $e) / 30.6001)
    $g1 = Int(30.6001 * $g)

    $tTimePlace.day = Int($c - $e - $g1)
    $tTimePlace.hour = ($jd - $jdi) * 24
    If $g <= 13 Then
        $tTimePlace.month = Int($g - 1)
    Else
        $tTimePlace.month = Int($g - 13)
    EndIf
    If $tTimePlace.month > 2 Then
        $tTimePlace.year = Int($d - 4716)
    Else
        $tTimePlace.year = Int($d - 4715)
    EndIf
EndFunc   ;==>JulianToDate

Func Julian($year, $month, $day)
    Local $a, $b, $c, $e
    If $month < 3 Then
        $year -= 1
        $month += 12
    EndIf
    If ($year > 1582 Or ($year = 1582 And $month > 10) Or ($year = 1582 And $month = 10 And $day > 15)) Then
        $a = ($year / 100)
        $b = Int(2 - $a + $a / 4)
    EndIf
    $c = Int(365.25 * $year)
    $e = Int(30.6001 * ($month + 1))
    Return $b + $c + $e + $day + 1720994.5
EndFunc   ;==>Julian

Func Sun_Position($j)
    Local $n, $x, $e, $l, $dl, $v, $m2, $i
    $n = 360 / 365.2422 * $j
    $i = Int($n / 360)
    $n = $n - $i * 360
    $x = $n - 3.762863
    If $x < 0 Then $x += 360
    $x *= $fRad
    $e = $x
    Do
        $dl = $e - 0.016718 * Sin($e) - $x
        $e = $e - $dl / (1 - 0.016718 * Cos($e))
    Until Abs($dl) < $fSmall_Float
    $v = 360 / $fPI * ATan(1.01686011182 * Tan($e / 2))
    $l = $v + 282.596403
    $i = Int($l / 360)
    $l = $l - $i * 360.0
    Return $l
EndFunc   ;==>Sun_Position

Func Moon_Position($j, $ls)
    Local $ms, $l, $mm, $n, $ev, $sms, $z, $x, $lm, $bm, $ae, $ec, $d, $ds, $as, $dm, $i
    $ms = 0.985647332099 * $j - 3.762863
    If ($ms < 0) Then $ms += 360.0
    $l = 13.176396 * $j + 64.975464
    $i = Int($l / 360)
    $l = $l - $i * 360.0
    If ($l < 0) Then $l += 360.0
    $mm = $l - 0.1114041 * $j - 349.383063
    $i = Int($mm / 360)
    $mm -= $i * 360.0
    $n = 151.950429 - 0.0529539 * $j
    $i = Int($n / 360)
    $n -= $i * 360.0
    $ev = 1.2739 * Sin((2 * ($l - $ls) - $mm) * $fRad)
    $sms = Sin($ms * $fRad)
    $ae = 0.1858 * $sms
    $mm += $ev - $ae - 0.37 * $sms
    $ec = 6.2886 * Sin($mm * $fRad)
    $l += $ev + $ec - $ae + 0.214 * Sin(2 * $mm * $fRad)
    $l = 0.6583 * Sin(2 * ($l - $ls) * $fRad) + $l
    Return $l
EndFunc   ;==>Moon_Position

Func Moon_Phase($year, $month, $day, $hour, ByRef $ip)
    Local $j = Julian($year, $month, $day + $hour / 24.0) - 2444238.5
    Local $ls = Sun_Position($j)
    Local $lm = Moon_Position($j, $ls)
    Local $t = $lm - $ls
    If ($t < 0) Then $t += 360
    $ip = BitAND(Int(($t + 22.5) / 45), 0x7)
    Return (1.0 - Cos(($lm - $ls) * $fRad)) / 2
EndFunc   ;==>Moon_Phase

Func NextDay(ByRef $y, ByRef $m, ByRef $d, $dd)
    Local $jd = Julian($y, $m, $d)
    $jd += $dd
    JulianToDate($tTimePlace, $jd)
    $y = $tTimePlace.year
    $m = $tTimePlace.month
    $d = $tTimePlace.day
EndFunc   ;==>NextDay

Func Main($y = @YEAR, $m = @MON)
    Local $d, $m0, $h, $i, $p, $ip = 0
    Local $step = 1, $begun = 0, $pmax = 0, $pmin = 1, $plast = 0
    Local $ymax, $mmax, $dmax, $hmax, $ymin, $mmin, $dmin, $hmin
    $d = 1
    $m0 = $m
    ConsoleWrite("tabulation of the phase of the moon for one month" & @CRLF & @CRLF)
    ConsoleWrite("Date       Time   Phase Segment" & @CRLF)
    Do
        For $h = 0 To 23 Step $step
            $p = Moon_Phase($y, $m, $d, $h, $ip)
            If $begun Then
                If ($p > $plast And $p > $pmax) Then
                    $pmax = $p
                    $ymax = ($y)
                    $mmax = ($m)
                    $dmax = ($d)
                    $hmax = ($h)
                ElseIf ($pmax) Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00          (fullest)\n", $ymax, $mmax, $dmax, $hmax))
                    $pmax = 0
                EndIf

                If $p < $plast And $p < $pmin Then
                    $pmin = $p
                    $ymin = ($y)
                    $mmin = ($m)
                    $dmin = ($d)
                    $hmin = ($h)
                ElseIf $pmin < 1 Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00          (newest)\n", $ymin, $mmin, $dmin, $hmin))
                    $pmin = 1.0
                EndIf
                If $h = 16 Then
                    ConsoleWrite(StringFormat("%04d/%02d/%02d %02d:00 %5.1f%%   (%d)\n", $y, $m, $d, $h, Floor($p * 1000 + 0.5) / 10, $ip))
                EndIf
            Else
                $begun = 1
            EndIf
            $plast = $p
        Next
        NextDay($y, $m, $d, 1)
        If ($m <> $m0) Then ExitLoop
    Until False
EndFunc   ;==>Main


Main()

 

Please don't send me any personal message and ask for support! I will not reply!

Selection of finest graphical examples at Codepen.io

The own fart smells best!
Her 'sikim hıyar' diyene bir avuç tuz alıp koşma!
¯\_(ツ)_/¯  ٩(●̮̮̃•̃)۶ ٩(-̮̮̃-̃)۶ૐ

Link to comment
Share on other sites

Thanks @LarsJ, you are right about the loop, I've just changed While - Wend with Do - Until and result given by the script are exactly the same as the compiled original c script. Fixed script in post #3
Thanks @UEZ for the nice overall rearrangement of the script, Interesting the use of the DllStructCreate

...I will modify this script to be a general purpose function

Thanks

Edited by Chimp

 

image.jpeg.9f1a974c98e9f77d824b358729b089b0.jpeg Chimp

small minds discuss people average minds discuss events great minds discuss ideas.... and use AutoIt....

Link to comment
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
 Share

  • Recently Browsing   0 members

    • No registered users viewing this page.
×
×
  • Create New...