1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using xWatchPhone;

namespace SonneMond
{
    public partial class SonneMond : Form
    {

        public SonneMond()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            SunMoonEngine sme = new SunMoonEngine(52,13,36,DateTime.Now);
            SunMoonEngine.SunMoonRiseSet smr = new SunMoonEngine.SunMoonRiseSet();
            label1.Text = Convert.ToString(smr.sunRise);
        }
    }
}					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:
228:
229:
230:
231:
232:
233:
234:
235:
236:
237:
238:
239:
240:
241:
242:
243:
244:
245:
246:
247:
248:
249:
250:
251:
252:
253:
254:
255:
256:
257:
258:
259:
260:
261:
262:
263:
264:
265:
266:
267:
268:
269:
270:
271:
272:
273:
274:
275:
276:
277:
278:
279:
280:
281:
282:
283:
284:
285:
286:
287:
288:
289:
290:
291:
292:
293:
294:
295:
296:
297:
298:
299:
300:
301:
302:
303:
304:
305:
306:
307:
308:
309:
310:
311:
312:
313:
314:
315:
316:
317:
318:
319:
320:
321:
322:
323:
324:
325:
326:
327:
328:
329:
330:
331:
332:
333:
334:
335:
336:
337:
338:
339:
340:
341:
342:
343:
344:
345:
346:
347:
348:
349:
350:
351:
352:
353:
354:
355:
356:
357:
358:
359:
360:
361:
362:
363:
364:
365:
366:
367:
368:
369:
370:
371:
372:
373:
374:
375:
376:
377:
378:
379:
380:
381:
382:
383:
384:
385:
386:
387:
388:
389:
390:
391:
392:
393:
394:
395:
396:
397:
398:
399:
400:
401:
402:
403:
404:
405:
406:
407:
408:
409:
410:
411:
412:
413:
414:
415:
416:
417:
418:
419:
420:
421:
422:
423:
424:
425:
426:
427:
428:
429:
430:
431:
432:
433:
434:
435:
436:
437:
438:
439:
440:
441:
442:
443:
444:
445:
446:
447:
448:
449:
450:
451:
452:
453:
454:
455:
456:
457:
458:
459:
460:
461:
462:
463:
464:
465:
466:
467:
468:
469:
470:
471:
472:
473:
474:
475:
476:
477:
478:
479:
480:
481:
482:
483:
484:
485:
486:
487:
488:
489:
490:
491:
492:
493:
494:
495:
496:
497:
498:
499:
500:
501:
502:
503:
504:
505:
506:
507:
508:
509:
510:
511:
512:
513:
514:
515:
516:
517:
518:
519:
520:
521:
522:
523:
524:
525:
526:
527:
528:
529:
530:
531:
532:
533:
534:
535:
536:
537:
538:
539:
540:
541:
542:
543:
544:
545:
546:
547:
548:
549:
550:
551:
552:
553:
554:
555:
556:
557:
558:
559:
560:
561:
562:
563:
564:
565:
566:
567:
568:
569:
570:
571:
572:
573:
574:
575:
576:
577:
578:
579:
580:
581:
582:
583:
584:
585:
586:
587:
588:
589:
590:
591:
592:
593:
594:
595:
596:
597:
598:
599:
600:
601:
602:
603:
604:
605:
606:
607:
608:
609:
610:
611:
612:
613:
614:
615:
616:
617:
618:
619:
620:
621:
622:
623:
624:
625:
626:
627:
628:
629:
630:
631:
632:
633:
634:
635:
636:
637:
638:
639:
640:
641:
642:
643:
644:
645:
646:
647:
648:
649:
650:
651:
652:
653:
654:
655:
656:
657:
658:
659:
660:
661:
662:
663:
664:
665:
666:
667:
668:
669:
670:
671:
672:
673:
674:
675:
676:
677:
678:
679:
680:
681:
682:
683:
684:
685:
686:
687:
688:
689:
690:
691:
692:
693:
694:
695:
696:
697:
698:
699:
700:
701:
702:
703:
704:
705:
706:
707:
708:
709:
710:
711:
712:
713:
714:
715:
716:
717:
718:
719:
720:
721:
722:
723:
724:
725:
726:
727:
728:
729:
730:
731:
732:
733:
734:
735:
736:
737:
738:
739:
740:
741:
742:
743:
744:
745:
746:
747:
748:
749:
750:
751:
752:
753:
754:
755:
756:
757:
758:
759:
760:
761:
762:
763:
764:
765:
766:
767:
768:
769:
770:
771:
772:
773:
774:
775:
776:
777:
778:
779:
780:
781:
782:
783:
784:
785:
786:
787:
788:
789:
790:
791:
792:
793:
794:
795:
796:
797:
798:
799:
800:
801:
802:
803:
804:
805:
806:
807:
808:
809:
810:
811:
812:
813:
814:
815:
816:
817:
818:
819:
820:
821:
822:
823:
824:
825:
826:
827:
828:
829:
830:
831:
832:
833:
834:
835:
836:
837:
838:
839:
840:
841:
842:
843:
844:
845:
846:
847:
848:
849:
850:
851:
852:
853:

/* 
  Entfernen Sie folgende Informationen auf keinen Fall: / Do not remove following text:
  Source code based on the javascript by Arnold Barmettler, www.astronomie.info / www.CalSky.com
  based on algorithms by Peter Duffett-Smith's great and easy book
  'Practical Astronomy with your Calculator'.
*/

using System;

namespace xWatchPhone
{
  public class SunMoonEngine
  {

    private struct coor
    {
      public double lon;
      public double lat;

      public double ra;
      public double dec;
      public double raGeocentric;
      public double decGeocentric;

      public double az;
      public double alt;

      public double x;
      public double y;
      public double z;        

      public double radius;
      public double diameter;
      public double distance;
      public double distanceTopocentric;
      public double decTopocentric;
      public double raTopocentric;

      public double anomalyMean;
      public double parallax;
      public double orbitLon;

      public double moonAge;
      public double phase;

      public string moonPhase;
      public string sign;
    }

    private struct riseset
    {
      public double transit;
      public double rise;
      public double set;

      public double cicilTwilightMorning;
      public double cicilTwilightEvening;

      public double nauticalTwilightMorning;
      public double nauticalTwilightEvening;

      public double astronomicalTwilightMorning;
      public double astronomicalTwilightEvening;
    }

    public struct SunMoonRiseSet
    {
      public DateTime DayOfValues;
      public double twilightMorning;
      public double sunRise;
      public double sunCulmination;
      public double sunSet;
      public double twilightEvening;

      public double moonRise;
      public double moonCulmination;
      public double moonSet;
      public double moonAge;     // 0 .. Neumond - 180 .. Vollmond - 359 .. Neumond
      public string moonPhase;
    }


    private double pi = 0;
    private double DEG = 0;
    private double RAD = 0;
    //private string empty = "--";

    private double lat;
    private double lon;
    private double seaLevel;
    private DateTime currDay;


    public SunMoonEngine()
    {
      pi = Math.PI;
      DEG = pi / 180.0;
      RAD = 180.0 / pi;

      lat = 52;
      lon = 13;
      seaLevel = 36;
      currDay = DateTime.Now;
    }

    public SunMoonEngine(double latitude, double longitude, double meterAboveSeaLevel, DateTime interestingDay)
    {
      pi = Math.PI;
      DEG = pi/180.0;
      RAD = 180.0/pi;

      lat = latitude;
      lon = longitude;
      seaLevel = meterAboveSeaLevel;
      currDay = interestingDay;
    }

    public SunMoonRiseSet InitEngine(double latitude, double longitude, double meterAboveSeaLevel, DateTime interestingDay)
    {
      lat = latitude;
      lon = longitude;
      seaLevel = meterAboveSeaLevel;
      currDay = interestingDay;

      return Compute();
    }


    private double sqr(double x) { return x * x; }

    // return integer value, closer to 0
    private int Int(double x) { return (x < 0) ? Convert.ToInt32(Math.Ceiling(x)) : Convert.ToInt32(Math.Floor(x)); }

    private double frac(double x) { return (x - Math.Floor(x)); }

    private double Mod(double a, double b) { return (a - Math.Floor(a / b) * b); }


    // Modulo PI
    private double Mod2Pi(double x) { return Mod(x, 2.0 * pi); }

    private double round100000(double x) { return Math.Round(100000.0 * x) / 100000.0; }
    private double round10000(double x) { return Math.Round(10000.0 * x) / 10000.0; }
    private double round1000(double x) { return Math.Round(1000.0 * x) / 1000.0; }
    private double round100(double x) { return Math.Round(100.0 * x) / 100.0; }
    private double round10(double x) { return Math.Round(10.0 * x) / 10.0; }



    private string HHMM(double hh)
    {
      if (hh == 0)
        return "";

      var m = frac(hh) * 60.0;
      var h = Int(hh);
      if (m >= 59.5) { h++; m -= 60.0; }
      m = Math.Round(m);

      var sH = "";

      if (h < 10) sH = "0";
      sH += h + ":";

      if (m < 10) sH += "0";
      sH += Convert.ToInt32(m);

      return sH;
    }

    //private string HHMMSS(double hh) 
    //{
    //    if (hh==0)
    //        return empty;
  
    //    var m = frac(hh)*60;
    //    var h = Int(hh);
    //    var s = frac(m)*60.0;
    //    var time="";
    //    m = Int(m);
    //    if (s>=59.5) { m++; s -=60.0; }
    //    if (m>=60)   { h++; m -=60; }
    //    s = Math.Round(s);
    //    if (h<10) time = "0"+h;
    //    else time=h.ToString();
    //    time += ":";
    //    if (m<10) time+="0"+m;
    //    else time+=m;
    //    time+=":";
    //    if (s<10) time+="0"+s;
    //    else time+=s;
    //    return time+" = "+round10000(hh);
    //}


    private string Sign(double lon)
    {
      string[] signs = {"Widder", "Stier", "Zwillinge", "Krebs", "Löwe", "Jungfrau", 
      "Waage", "Skorpion", "Schütze", "Steinbock", "Wassermann", "Fische"};

      return (signs[Convert.ToInt16(Math.Floor(lon * RAD / 30))]);
    }


    // Calculate Julian date: valid only from 1.3.1901 to 28.2.2100
    private double CalcJD(int day, int month, int year)
    {
      var jd = 2415020.5 - 64; // 1.1.1900 - correction of algorithm
      if (month <= 2) { year--; month += 12; }
      jd += Int((year - 1900) * 365.25);
      jd += Int(30.6001 * (1 + month));

      return jd + day;
    }


    // Julian Date to Greenwich Mean Sidereal Time
    private double GMST(double JD)
    {
      var UT = frac(JD - 0.5) * 24.0; // UT in hours
      JD = Math.Floor(JD - 0.5) + 0.5;   // JD at 0 hours UT
      var T = (JD - 2451545.0) / 36525.0;
      var T0 = 6.697374558 + T * (2400.051336 + T * 0.000025862);

      return Mod(T0 + UT * 1.002737909, 24.0);
    }


    // Convert Greenweek mean sidereal time to UT
    private double GMST2UT(double JD, double gmst)
    {
      JD = Math.Floor(JD - 0.5) + 0.5;   // JD at 0 hours UT
      var T = (JD - 2451545.0) / 36525.0;
      var T0 = Mod(6.697374558 + T * (2400.051336 + T * 0.000025862), 24.0);
      //var UT = 0.9972695663*Mod((gmst-T0), 24.);
      var UT = 0.9972695663 * ((gmst - T0));

      return UT;
    }


    // Local Mean Sidereal Time, geographical longitude in radians, East is positive
    private double GMST2LMST(double gmst, double lon)
    {
      var lmst = Mod(gmst + RAD * lon / 15, 24.0);
      return lmst;
    }


    // Transform ecliptical coordinates (lon/lat) to equatorial coordinates (RA/dec)
    private coor Ecl2Equ(coor c0, double TDT)
    {
      var T = (TDT - 2451545.0) / 36525.0; // Epoch 2000 January 1.5
      var eps = (23.0 + (26 + 21.45 / 60.0) / 60.0 + T * (-46.815 + T * (-0.0006 + T * 0.00181)) / 3600.0) * DEG;
      var coseps = Math.Cos(eps);
      var sineps = Math.Sin(eps);

      var sinlon = Math.Sin(c0.lon);
      coor c1 = c0;
      c1.ra = Mod2Pi(Math.Atan2((sinlon * coseps - Math.Tan(c0.lat) * sineps), Math.Cos(c0.lon)));
      c1.dec = Math.Asin(Math.Sin(c0.lat) * coseps + Math.Cos(c0.lat) * sineps * sinlon);

      return c1;
    }


    // Transform equatorial coordinates (RA/Dec) to horizonal coordinates (azimuth/altitude)
    // Refraction is ignored
    private coor Equ2Altaz(coor c0, double TDT, double geolat, double lmst)
    {
      var cosdec = Math.Cos(c0.dec);
      var sindec = Math.Sin(c0.dec);
      var lha = lmst - c0.ra;
      var coslha = Math.Cos(lha);
      var sinlha = Math.Sin(lha);
      var coslat = Math.Cos(geolat);
      var sinlat = Math.Sin(geolat);

      var N = -cosdec * sinlha;
      var D = sindec * coslat - cosdec * coslha * sinlat;
      coor c1 = c0;
      c1.az = Mod2Pi(Math.Atan2(N, D));
      c1.alt = Math.Asin(sindec * sinlat + cosdec * coslha * coslat);

      return c1;
    }


    // Transform geocentric equatorial coordinates (RA/Dec) to topocentric equatorial coordinates
    private coor GeoEqu2TopoEqu(coor c0, coor observer, double lmst)
    {
      var cosdec = Math.Cos(c0.dec);
      var sindec = Math.Sin(c0.dec);
      var coslst = Math.Cos(lmst);
      var sinlst = Math.Sin(lmst);
      var coslat = Math.Cos(observer.lat); // we should use geocentric latitude, not geodetic latitude
      var sinlat = Math.Sin(observer.lat);
      var rho = observer.radius; // observer-geocenter in Kilometer

      var x = c0.distance * cosdec * Math.Cos(c0.ra) - rho * coslat * coslst;
      var y = c0.distance * cosdec * Math.Sin(c0.ra) - rho * coslat * sinlst;
      var z = c0.distance * sindec - rho * sinlat;

      var c1 = c0;
      c1.distanceTopocentric = Math.Sqrt(x * x + y * y + z * z);
      c1.decTopocentric = Math.Asin(z / c1.distanceTopocentric);
      c1.raTopocentric = Mod2Pi(Math.Atan2(y, x));

      return c1;
    }


    // Calculate cartesian from polar coordinates
    private coor EquPolar2Cart(double lon, double lat, double distance)
    {
      coor c1 = new coor();
      var rcd = Math.Cos(lat) * distance;
      c1.x = rcd * Math.Cos(lon);
      c1.y = rcd * Math.Sin(lon);
      c1.z = distance * Math.Sin(lat);
      return c1;
    }


    // Calculate observers cartesian equatorial coordinates (x,y,z in celestial frame) 
    // from geodetic coordinates (longitude, latitude, height above WGS84 ellipsoid)
    // Currently only used to calculate distance of a body from the observer
    private coor Observer2EquCart(double lon, double lat, double height, double gmst)
    {
      var flat = 298.257223563;        // WGS84 flatening of earth
      var aearth = 6378.137;           // GRS80/WGS84 semi major axis of earth ellipsoid
      coor c1 = new coor();
      // Calculate geocentric latitude from geodetic latitude
      var co = Math.Cos(lat);
      var si = Math.Sin(lat);
      var fl = 1.0 - 1.0 / flat;
      fl = fl * fl;
      si = si * si;
      var u = 1.0 / Math.Sqrt(co * co + fl * si);
      var a = aearth * u + height;
      var b = aearth * fl * u + height;
      var radius = Math.Sqrt(a * a * co * co + b * b * si); // geocentric distance from earth center
      c1.y = Math.Acos(a * co / radius); // geocentric latitude, rad
      c1.x = lon; // longitude stays the same
      if (lat < 0.0) { c1.y = -c1.y; } // adjust sign
      c1 = EquPolar2Cart(c1.x, c1.y, radius); // convert from geocentric polar to geocentric cartesian, with regard to Greenwich
      // rotate around earth's polar axis to align coordinate system from Greenwich to vernal equinox
      var x = c1.x;
      var y = c1.y;
      var rotangle = gmst / 24 * 2 * pi; // sideral time gmst given in hours. Convert to radians
      c1.x = x * Math.Cos(rotangle) - y * Math.Sin(rotangle);
      c1.y = x * Math.Sin(rotangle) + y * Math.Cos(rotangle);
      c1.radius = radius;
      c1.lon = lon;
      c1.lat = lat;
      return c1;
    }


    // Calculate coordinates for Sun
    // Coordinates are accurate to about 10s (right ascension) 
    // and a few minutes of arc (declination)
    private coor SunPosition(double TDT, double geolat, double lmst)
    {
      var D = TDT - 2447891.5;

      var eg = 279.403303 * DEG;
      var wg = 282.768422 * DEG;
      var e = 0.016713;
      var a = 149598500; // km
      var diameter0 = 0.533128 * DEG; // angular diameter of Moon at a distance

      var MSun = 360 * DEG / 365.242191 * D + eg - wg;
      var nu = MSun + 360.0 * DEG / pi * e * Math.Sin(MSun);

      var sunCoor = new coor();
      sunCoor.lon = Mod2Pi(nu + wg);
      sunCoor.lat = 0;
      sunCoor.anomalyMean = MSun;

      sunCoor.distance = (1 - sqr(e)) / (1 + e * Math.Cos(nu)); // distance in astronomical units
      sunCoor.diameter = diameter0 / sunCoor.distance; // angular diameter in radians
      sunCoor.distance *= a;                         // distance in km
      sunCoor.parallax = 6378.137 / sunCoor.distance;  // horizonal parallax

      var sunCoor1 = Ecl2Equ(sunCoor, TDT);

      // Calculate horizonal coordinates of sun
      var sunCoor2 = Equ2Altaz(sunCoor1, TDT, geolat, lmst);
      sunCoor2.sign = Sign(sunCoor.lon);

      return sunCoor2;
    }


    // Calculate data and coordinates for the Moon
    // Coordinates are accurate to about 1/5 degree (in ecliptic coordinates)
    private coor MoonPosition(coor sunCoor, double TDT, coor observer, double lmst)
    {
      var D = TDT - 2447891.5;

      // Mean Moon orbit elements as of 1990.0
      var l0 = 318.351648 * DEG;
      var P0 = 36.340410 * DEG;
      var N0 = 318.510107 * DEG;
      var i = 5.145396 * DEG;
      var e = 0.054900;
      var a = 384401;                // km
      var diameter0 = 0.5181 * DEG;     // angular diameter of Moon at a distance
      var parallax0 = 0.9507 * DEG;     // parallax at distance a

      var l = 13.1763966 * DEG * D + l0;
      var MMoon = l - 0.1114041 * DEG * D - P0; // Moon's mean anomaly M
      var N = N0 - 0.0529539 * DEG * D;       // Moon's mean ascending node longitude
      var C = l - sunCoor.lon;
      var Ev = 1.2739 * DEG * Math.Sin(2 * C - MMoon);
      var Ae = 0.1858 * DEG * Math.Sin(sunCoor.anomalyMean);
      var A3 = 0.37 * DEG * Math.Sin(sunCoor.anomalyMean);
      var MMoon2 = MMoon + Ev - Ae - A3;            // corrected Moon anomaly
      var Ec = 6.2886 * DEG * Math.Sin(MMoon2);   // equation of centre
      var A4 = 0.214 * DEG * Math.Sin(2 * MMoon2);
      var l2 = l + Ev + Ec - Ae + A4;                 // corrected Moon's longitude
      var V = 0.6583 * DEG * Math.Sin(2 * (l2 - sunCoor.lon));
      var l3 = l2 + V;                          // true orbital longitude;

      var N2 = N - 0.16 * DEG * Math.Sin(sunCoor.anomalyMean);

      var moonCoor = new coor();
      moonCoor.lon = Mod2Pi(N2 + Math.Atan2(Math.Sin(l3 - N2) * Math.Cos(i), Math.Cos(l3 - N2)));
      moonCoor.lat = Math.Asin(Math.Sin(l3 - N2) * Math.Sin(i));
      moonCoor.orbitLon = l3;

      moonCoor = Ecl2Equ(moonCoor, TDT);
      // relative distance to semi mayor axis of lunar oribt
      moonCoor.distance = (1 - sqr(e)) / (1 + e * Math.Cos(MMoon2 + Ec));
      moonCoor.diameter = diameter0 / moonCoor.distance; // angular diameter in radians
      moonCoor.parallax = parallax0 / moonCoor.distance; // horizontal parallax in radians
      moonCoor.distance *= a; // distance in km

      // Calculate horizonal coordinates of sun
      // transform geocentric coordinates into topocentric (==observer based) coordinates
      moonCoor = GeoEqu2TopoEqu(moonCoor, observer, lmst);
      moonCoor.raGeocentric = moonCoor.ra; // backup geocentric coordinates
      moonCoor.decGeocentric = moonCoor.dec;
      moonCoor.ra = moonCoor.raTopocentric;
      moonCoor.dec = moonCoor.decTopocentric;
      moonCoor = Equ2Altaz(moonCoor, TDT, observer.lat, lmst); // now ra and dec are topocentric

      // Age of Moon in radians since New Moon (0) - Full Moon (pi)
      moonCoor.moonAge = Mod2Pi(l3 - sunCoor.lon);
      moonCoor.phase = 0.5 * (1 - Math.Cos(moonCoor.moonAge)); // Moon phase, 0-1

      //string[] phases = { "Neumond", "Zunehmende Sichel", "Erstes Viertel", "Zunehmender Mond", "Vollmond", "Abnehmender Mond", "Letztes Viertel", "Abnehmende Sichel", "Neumond" };
      string[] phases = { "New moon", "Increasing crescent", "1st quarter", "Increasing Moon", "Full moon", "Decreasing Moon", "Last quarter", "decreasing crescent", "New moon" };
      var mainPhase = 1.0 / 29.53 * 360 * DEG; // show 'Newmoon, 'Quarter' for +/-1 day arond the actual event
      var p = Mod(moonCoor.moonAge, 90.0 * DEG);

      if (p < mainPhase || p > 90 * DEG - mainPhase)
        p = 2 * Math.Round(moonCoor.moonAge / (90.0 * DEG));
      else
        p = 2 * Math.Floor(moonCoor.moonAge / (90.0 * DEG)) + 1;

      moonCoor.moonPhase = phases[Convert.ToInt16(p)];

      moonCoor.sign = Sign(moonCoor.lon);

      return moonCoor;
    }


    // Rough refraction formula using standard atmosphere: 1015 mbar and 10°C
    // Input true altitude in radians, Output: increase in altitude in degrees
    private double Refraction(double alt)
    {
      var pressure = 1015;
      var temperature = 10;
      var altdeg = alt * RAD;

      if (altdeg < -2 || altdeg >= 90)
        return 0;

      if (altdeg > 15)
        return 0.00452 * pressure / ((273 + temperature) * Math.Tan(alt));

      var y = alt;
      var D = 0.0;
      var P = (pressure - 80.0) / 930.0;
      var Q = 0.0048 * (temperature - 10.0);
      var y0 = y;
      var D0 = D;
      var N = 0.0;

      for (var i = 0; i < 3; i++)
      {
        N = y + (7.31 / (y + 4.4));
        N = 1.0 / Math.Tan(N * DEG);
        D = N * P / (60.0 + Q * (N + 39.0));
        N = y - y0;
        y0 = D - D0 - N;
        N = ((N != 0.0) && (y0 != 0.0)) ? y - N * (alt + D - y) / y0 : alt + D;
        y0 = y;
        D0 = D;
        y = N;
      }
      return D; // Hebung durch Refraktion in radians
    }


    // returns Greenwich sidereal time (hours) of time of rise 
    // and set of object with coordinates coor.ra/coor.dec
    // at geographic position lon/lat (all values in radians)
    // Correction for refraction and semi-diameter/parallax of body is taken care of in function RiseSet
    // h is used to calculate the twilights. It gives the required elevation of the disk center of the sun
    private riseset GMSTRiseSet(coor c0, double lon, double lat, double h)
    {
      //var h = (h == null) ? 0.0 : h; // set default value
      // var tagbogen = Math.acos(-Math.tan(lat)*Math.tan(coor.dec)); // simple formula if twilight is not required
      var tagbogen = Math.Acos((Math.Sin(h) - Math.Sin(lat) * Math.Sin(c0.dec)) / (Math.Cos(lat) * Math.Cos(c0.dec)));

      riseset r1 = new SunMoonEngine.riseset();
      r1.transit = RAD / 15 * (+c0.ra - lon);
      r1.rise = 24.0 + RAD / 15 * (-tagbogen + c0.ra - lon); // calculate GMST of rise of object
      r1.set = RAD / 15 * (+tagbogen + c0.ra - lon); // calculate GMST of set of object

      // using the modulo function Mod, the day number goes missing. This may get a problem for the moon
      r1.transit = Mod(r1.transit, 24);
      r1.rise = Mod(r1.rise, 24);
      r1.set = Mod(r1.set, 24);

      return r1;
    }


    // Find GMST of rise/set of object from the two calculates 
    // (start)points (day 1 and 2) and at midnight UT(0)
    private double InterpolateGMST(double gmst0, double gmst1, double gmst2, double timefactor)
    {
      return (timefactor * 24.07 * gmst1 - gmst0 * (gmst2 - gmst1)) / (timefactor * 24.07 + gmst1 - gmst2);
    }


    // JD is the Julian Date of 0h UTC time (midnight)
    private riseset RiseSet(double jd0UT, coor c1, coor c2, double lon, double lat, double timeinterval, double altitude)
    {
      // altitude of sun center: semi-diameter, horizontal parallax and (standard) refraction of 34'
      var alt = 0.0; // calculate 
      //var altitude = (altitude == null) ? 0. : altitude; // set default value

      // true height of sun center for sunrise and set calculation. Is kept 0 for twilight (ie. altitude given):
      if (altitude == 0) alt = 0.5 * c1.diameter - c1.parallax + 34.0 / 60 * DEG;

      var rise1 = GMSTRiseSet(c1, lon, lat, altitude);
      var rise2 = GMSTRiseSet(c2, lon, lat, altitude);

      var rise = new riseset();

      // unwrap GMST in case we move across 24h -> 0h
      if (rise1.transit > rise2.transit && Math.Abs(rise1.transit - rise2.transit) > 18) rise2.transit += 24;
      if (rise1.rise > rise2.rise && Math.Abs(rise1.rise - rise2.rise) > 18) rise2.rise += 24;
      if (rise1.set > rise2.set && Math.Abs(rise1.set - rise2.set) > 18) rise2.set += 24;
      var T0 = GMST(jd0UT);
      //  var T02 = T0-zone*1.002738; // Greenwich sidereal time at 0h time zone (zone: hours)

      // Greenwich sidereal time for 0h at selected longitude
      var T02 = T0 - lon * RAD / 15 * 1.002738; if (T02 < 0) T02 += 24;

      if (rise1.transit < T02) { rise1.transit += 24; rise2.transit += 24; }
      if (rise1.rise < T02) { rise1.rise += 24; rise2.rise += 24; }
      if (rise1.set < T02) { rise1.set += 24; rise2.set += 24; }

      // Refraction and Parallax correction
      var decMean = 0.5 * (c1.dec + c2.dec);
      var psi = Math.Acos(Math.Sin(lat) / Math.Cos(decMean));
      var y = Math.Asin(Math.Sin(alt) / Math.Sin(psi));
      var dt = 240 * RAD * y / Math.Cos(decMean) / 3600; // time correction due to refraction, parallax

      rise.transit = GMST2UT(jd0UT, InterpolateGMST(T0, rise1.transit, rise2.transit, timeinterval));
      rise.rise = GMST2UT(jd0UT, InterpolateGMST(T0, rise1.rise, rise2.rise, timeinterval) - dt);
      rise.set = GMST2UT(jd0UT, InterpolateGMST(T0, rise1.set, rise2.set, timeinterval) + dt);

      return rise;
    }


    // Find (local) time of sunrise and sunset, and twilights
    // JD is the Julian Date of 0h local time (midnight)
    // Accurate to about 1-2 minutes
    // recursive: 1 - calculate rise/set in UTC in a second run
    // recursive: 0 - find rise/set on the current local day. This is set when doing the first call to this function
    private riseset SunRise(double JD, double deltaT, double lon, double lat, double zone, bool recursive)
    {
      var jd0UT = Math.Floor(JD - 0.5) + 0.5;   // JD at 0 hours UT
      var coor1 = SunPosition(jd0UT + deltaT / 24.0 / 3600.0, 0, 0);
      var coor2 = SunPosition(jd0UT + 1.0 + deltaT / 24.0 / 3600.0, 0, 0); // calculations for next day's UTC midnight

      var risetemp = new riseset();
      var rise = new riseset();
      // rise/set time in UTC. 
      rise = RiseSet(jd0UT, coor1, coor2, lon, lat, 1, 0);
      if (!recursive)
      { // check and adjust to have rise/set time on local calendar day
        if (zone > 0)
        {
          // rise time was yesterday local time -> calculate rise time for next UTC day
          if (rise.rise >= 24 - zone || rise.transit >= 24 - zone || rise.set >= 24 - zone)
          {
            risetemp = SunRise(JD + 1, deltaT, lon, lat, zone, true);
            if (rise.rise >= 24 - zone) rise.rise = risetemp.rise;
            if (rise.transit >= 24 - zone) rise.transit = risetemp.transit;
            if (rise.set >= 24 - zone) rise.set = risetemp.set;
          }
        }
        else if (zone < 0)
        {
          // rise time was yesterday local time -> calculate rise time for next UTC day
          if (rise.rise < -zone || rise.transit < -zone || rise.set < -zone)
          {
            risetemp = SunRise(JD - 1, deltaT, lon, lat, zone, true);
            if (rise.rise < -zone) rise.rise = risetemp.rise;
            if (rise.transit < -zone) rise.transit = risetemp.transit;
            if (rise.set < -zone) rise.set = risetemp.set;
          }
        }

        rise.transit = Mod(rise.transit + zone, 24.0);
        rise.rise = Mod(rise.rise + zone, 24.0);
        rise.set = Mod(rise.set + zone, 24.0);

        // Twilight calculation
        // civil twilight time in UTC. 
        risetemp = RiseSet(jd0UT, coor1, coor2, lon, lat, 1, -6.0 * DEG);
        rise.cicilTwilightMorning = Mod(risetemp.rise + zone, 24.0);
        rise.cicilTwilightEvening = Mod(risetemp.set + zone, 24.0);

        // nautical twilight time in UTC. 
        risetemp = RiseSet(jd0UT, coor1, coor2, lon, lat, 1, -12.0 * DEG);
        rise.nauticalTwilightMorning = Mod(risetemp.rise + zone, 24.0);
        rise.nauticalTwilightEvening = Mod(risetemp.set + zone, 24.0);

        // astronomical twilight time in UTC. 
        risetemp = RiseSet(jd0UT, coor1, coor2, lon, lat, 1, -18.0 * DEG);
        rise.astronomicalTwilightMorning = Mod(risetemp.rise + zone, 24.0);
        rise.astronomicalTwilightEvening = Mod(risetemp.set + zone, 24.0);
      }
      return rise;
    }



    // Find local time of moonrise and moonset
    // JD is the Julian Date of 0h local time (midnight)
    // Accurate to about 5 minutes or better
    // recursive: 1 - calculate rise/set in UTC
    // recursive: 0 - find rise/set on the current local day (set could also be first)
    // returns 0.000000000 for moonrise/set does not occur on selected day
    private riseset MoonRise(double JD, double deltaT, double lon, double lat, double zone, bool recursive)
    {
      var timeinterval = 0.5;

      var jd0UT = Math.Floor(JD - 0.5) + 0.5;   // JD at 0 hours UT
      var suncoor1 = SunPosition(jd0UT + deltaT / 24.0 / 3600.0, 0, 0);
      var coor1 = MoonPosition(suncoor1, jd0UT + deltaT / 24.0 / 3600.0, new coor(), 0);

      var suncoor2 = SunPosition(jd0UT + timeinterval + deltaT / 24.0 / 3600.0, 0, 0); // calculations for noon
      // calculations for next day's midnight
      var coor2 = MoonPosition(suncoor2, jd0UT + timeinterval + deltaT / 24.0 / 3600.0, new coor(), 0);

      var risetemp = new riseset();
      var rise = new riseset();

      // rise/set time in UTC, time zone corrected later.
      // Taking into account refraction, semi-diameter and parallax
      rise = RiseSet(jd0UT, coor1, coor2, lon, lat, timeinterval, 0);

      if (!recursive)
      { // check and adjust to have rise/set time on local calendar day
        if (zone > 0)
        {
          // recursive call to MoonRise returns events in UTC
          var riseprev = MoonRise(JD - 1.0, deltaT, lon, lat, zone, true);

          // recursive call to MoonRise returns events in UTC
          //risenext = MoonRise(JD+1, deltaT, lon, lat, zone, 1);
          //alert("yesterday="+riseprev.transit+"  today="+rise.transit+" tomorrow="+risenext.transit);
          //alert("yesterday="+riseprev.rise+"  today="+rise.rise+" tomorrow="+risenext.rise);
          //alert("yesterday="+riseprev.set+"  today="+rise.set+" tomorrow="+risenext.set);

          if (rise.transit >= 24.0 - zone || rise.transit < -zone)
          { // transit time is tomorrow local time
            if (riseprev.transit < 24.0 - zone) rise.transit = 0.000000000; // there is no moontransit today
            else rise.transit = riseprev.transit;
          }

          if (rise.rise >= 24.0 - zone || rise.rise < -zone)
          { // transit time is tomorrow local time
            if (riseprev.rise < 24.0 - zone) rise.rise = 0.000000000; // there is no moontransit today
            else rise.rise = riseprev.rise;
          }

          if (rise.set >= 24.0 - zone || rise.set < -zone)
          { // transit time is tomorrow local time
            if (riseprev.set < 24.0 - zone) rise.set = 0.000000000; // there is no moontransit today
            else rise.set = riseprev.set;
          }

        }
        else if (zone < 0)
        {
          // rise/set time was tomorrow local time -> calculate rise time for former UTC day
          if (rise.rise < -zone || rise.set < -zone || rise.transit < -zone)
          {
            risetemp = MoonRise(JD + 1.0, deltaT, lon, lat, zone, true);

            if (rise.rise < -zone)
            {
              if (risetemp.rise > -zone) rise.rise = 0.000000000; // there is no moonrise today
              else rise.rise = risetemp.rise;
            }

            if (rise.transit < -zone)
            {
              if (risetemp.transit > -zone) rise.transit = 0.000000000; // there is no moonset today
              else rise.transit = risetemp.transit;
            }

            if (rise.set < -zone)
            {
              if (risetemp.set > -zone) rise.set = 0.000000000; // there is no moonset today
              else rise.set = risetemp.set;
            }

          }
        }

        if (rise.rise != 0.000000000) rise.rise = Mod(rise.rise + zone, 24.0);          // correct for time zone, if time is valid
        if (rise.transit != 0.000000000) rise.transit = Mod(rise.transit + zone, 24.0); // correct for time zone, if time is valid
        if (rise.set != 0.000000000) rise.set = Mod(rise.set + zone, 24.0);             // correct for time zone, if time is valid
      }
      return rise;
    }


    public double GetCurrLMST()
    {
      var JD0 = CalcJD(DateTime.Now.Day, DateTime.Now.Month, DateTime.Now.Year);
      var zone = (DateTime.Now - DateTime.UtcNow).Hours;
      var JD = JD0
        + (DateTime.Now.Hour - zone + DateTime.Now.Minute / 60.0
        + DateTime.Now.Second / 3600.0) / 24.0;

      var longitude = lon * DEG;      // longitude of observer

      var gmst = GMST(JD);
      return GMST2LMST(gmst, longitude);
    }

    public SunMoonRiseSet Compute()
    {

      if (currDay.Year <= 1900 || currDay.Year >= 2100)
        throw new Exception("Diese Berechnung erlaubt nur eine Zeitperiode 1901-2099!");

      var calculation = new SunMoonRiseSet();

            calculation.DayOfValues = currDay;

      var JD0 = CalcJD(currDay.Day, currDay.Month, currDay.Year);
      var zone = (DateTime.Now - DateTime.UtcNow).Hours;
      var JD = JD0
        + (currDay.Hour - zone + currDay.Minute / 60.0
        + currDay.Second / 3600.0) / 24.0;
      var DeltaT = 65.0;              // deltaT - difference among 'earth center' versus 'observered' time (TDT-UT), in seconds
      var TDT = JD + DeltaT / 24.0 / 3600.0;

      var latitude = lat * DEG;       // geodetic latitude of observer on WGS84
      var longitude = lon * DEG;      // longitude of observer
      var height = seaLevel * 0.001;  // altiude of observer in meters above WGS84 ellipsoid (and converted to kilometers)

      var gmst = GMST(JD);
      var lmst = GMST2LMST(gmst, longitude);

      var sunCoor = SunPosition(TDT, latitude, lmst * 15.0 * DEG);                // Calculate data for the Sun at given time

      //form.JD.value = round100000(JD);
      //form.GMST.value = HHMMSS(gmst);
      //form.LMST.value = HHMMSS(lmst);

      //if (eval(form.Minute.value)<10) form.Minute.value = "0"+eval(form.Minute.value);
      //if (eval(form.Month.value)<10) form.Month.value = "0"+eval(form.Month.value);

      //var SunLon = round1000(sunCoor.lon * RAD);
      //form.SunRA.value = HHMM(sunCoor.ra * RAD / 15);
      //var SunDec = round1000(sunCoor.dec * RAD);
      //var SunAz = round100(sunCoor.az * RAD);
      //var SunAlt = round10(sunCoor.alt * RAD + Refraction(sunCoor.alt));  // including refraction

      //form.SunSign.value = sunCoor.sign;
      //form.SunDiameter.value = round100(sunCoor.diameter*RAD*60.); // angular diameter in arc seconds
      //form.SunDistance.value = round10(sunCoor.distance);

      // Calculate distance from the observer (on the surface of earth) to the center of the sun
      //var sunCart = EquPolar2Cart(sunCoor.ra, sunCoor.dec, sunCoor.distance);
      //form.SunDistanceObserver.value = round10( Math.Sqrt( sqr(sunCart.x-observerCart.x) + sqr(sunCart.y-observerCart.y) + sqr(sunCart.z-observerCart.z) ));

      // JD0: JD of 0h UTC time
      var sunRise = SunRise(JD0, DeltaT, longitude, latitude, zone, false);

      calculation.sunCulmination = sunRise.transit;
      calculation.sunRise = sunRise.rise;
      calculation.sunSet = sunRise.set;

      calculation.twilightMorning = sunRise.cicilTwilightMorning;
      calculation.twilightEvening = sunRise.cicilTwilightEvening;
      //form.SunNauticalTwilightMorning.value    = HHMM(sunRise.nauticalTwilightMorning);
      //form.SunNauticalTwilightEvening.value    = HHMM(sunRise.nauticalTwilightEvening);
      //form.SunAstronomicalTwilightMorning.value    = HHMM(sunRise.astronomicalTwilightMorning);
      //form.SunAstronomicalTwilightEvening.value    = HHMM(sunRise.astronomicalTwilightEvening);



      var observerCart = Observer2EquCart(longitude, latitude, height, gmst);     // geocentric cartesian coordinates of observer
      var moonCoor = MoonPosition(sunCoor, TDT, observerCart, lmst * 15.0 * DEG); // Calculate data for the Moon at given time

      var MoonLon = round1000(moonCoor.lon*RAD);
      var MoonLat = round1000(moonCoor.lat*RAD);
      var MoonRA = HHMM(moonCoor.ra*RAD/15.0);
      var MoonDec = round1000(moonCoor.dec*RAD);
      var MoonAz = round100(moonCoor.az*RAD);
      var MoonAlt = round10(moonCoor.alt*RAD+Refraction(moonCoor.alt));  // including refraction

      calculation.moonAge = round1000(moonCoor.moonAge * RAD);

      var MoonPhaseNumber = round1000(moonCoor.phase);
      calculation.moonPhase = moonCoor.moonPhase;

      //form.MoonSign.value     = moonCoor.sign;
      var MoonDistance = round10(moonCoor.distance);
      var MoonDiameter = round100(moonCoor.diameter*RAD*60.0); // angular diameter in arc seconds

      // Calculate distance from the observer (on the surface of earth) to the center of the moon
      //var moonCart = EquPolar2Cart(moonCoor.raGeocentric, moonCoor.decGeocentric, moonCoor.distance);
      //form.MoonDistanceObserver.value = round10( Math.Sqrt( sqr(moonCart.x-observerCart.x) + sqr(moonCart.y-observerCart.y) + sqr(moonCart.z-observerCart.z) ));

      var moonRise = MoonRise(JD0, DeltaT, longitude, latitude, zone, false);

      calculation.moonCulmination = moonRise.transit;
      calculation.moonRise = moonRise.rise;
      calculation.moonSet = moonRise.set;

      return calculation;
    }
  }
}					

1:
2:

SunMoonEngine.SunMoonRiseSet smr = new SunMoonEngine.SunMoonRiseSet();
label1.Text = Convert.ToString(smr.sunRise);					

1:
2:
3:

SunMoonEngine sme = new SunMoonEngine(52,13,36,DateTime.Now);
SunMoonRiseSet smr = sme.Compute();
label1.Text = smr.sunRise.ToString(); // finde ich persönlich besser als Convert.ToString(...)					

1:
2:

int hour = (int)smr.sunRise;
int minute = (int)((smr.sunRise - (double)hour) * 60.0 + 0.5);