program generate_moon use, intrinsic :: iso_fortran_env, only: int32, real64 implicit none ! Parameter für das Bild (1000 x 750 Pixel) integer(int32), parameter :: width = 1000 integer(int32), parameter :: height = 750 real(real64), parameter :: pi = 3.14159265358979323846_real64 ! Standorte character(len=32) :: loc_name(2) real(real64) :: lat(2), lon(2) integer(int32) :: i loc_name(1) = "MINGA" lat(1) = 48.173333_real64 lon(1) = 11.635000_real64 loc_name(2) = "GAGGERS" lat(2) = 48.333333_real64 lon(2) = 11.250000_real64 do i = 1, 2 call render_moon_image(trim(loc_name(i)), lat(i), lon(i)) end do contains subroutine render_moon_image(location_name, latitude, longitude) character(len=*), intent(in) :: location_name real(real64), intent(in) :: latitude, longitude character(len=64) :: filename_ppm, filename_raw_png, filename_final_png character(len=4096) :: sys_cmd character(len=8) :: dt_date, dt_time character(len=5) :: dt_zone character(len=2048) :: overlay_text character(len=128) :: line_loc, line_zeit, line_phase, line_entf, line_az, line_alt integer(int32) :: img(width, height, 3) integer(int32) :: x, y, cx, cy, radius, fov_radius real(real64) :: dx, dy, dist, noise, brightness ! Datum, Zeit & Zeitzone integer :: year, month, day, hour, minute, second integer :: z_h, z_m real(real64) :: utc_offset_hours, julian_day, alt, az, phase_pct, dist_km character(len=32) :: phase_str ! Dateinamen festlegen filename_final_png = lowercase(location_name) // ".png" filename_raw_png = lowercase(location_name) // "_raw.png" filename_ppm = lowercase(location_name) // ".ppm" ! Systemzeit & Zeitzonen-Offset abfragen call date_and_time(date=dt_date, time=dt_time, zone=dt_zone) read(dt_date(1:4), '(I4)') year read(dt_date(5:6), '(I2)') month read(dt_date(7:8), '(I2)') day read(dt_time(1:2), '(I2)') hour read(dt_time(3:4), '(I2)') minute read(dt_time(5:6), '(I2)') second ! Zeitzonen-Offset in Stunden berechnen (z.B. +0200 -> +2.0) read(dt_zone(2:3), '(I2)') z_h read(dt_zone(4:5), '(I2)') z_m if (dt_zone(1:1) == '-') then utc_offset_hours = -1.0_real64 * (real(z_h, real64) + real(z_m, real64)/60.0_real64) else utc_offset_hours = 1.0_real64 * (real(z_h, real64) + real(z_m, real64)/60.0_real64) end if ! Astronomische Daten mit UTC berechnen julian_day = get_julian_day_utc(year, month, day, hour, minute, second, utc_offset_hours) call calculate_moon_position(julian_day, latitude, longitude, alt, az, phase_pct, dist_km, phase_str) ! --- GRAFIK RENDERN --- img = 0 cx = width / 2 cy = height / 2 + 30 fov_radius = 280 radius = 160 ! 8x24 Fernglas-Vergrößerung do y = 1, height do x = 1, width dx = real(x - cx, real64) dy = real(y - cy, real64) dist = sqrt(dx*dx + dy*dy) if (dist <= real(fov_radius, real64)) then img(x, y, 1) = 12 img(x, y, 2) = 15 img(x, y, 3) = 22 if (dist <= real(radius, real64)) then noise = sin(dx*0.15_real64)*cos(dy*0.15_real64) + sin(dx*0.05_real64 + dy*0.08_real64) brightness = 180.0_real64 + noise * 50.0_real64 brightness = brightness * sqrt(1.0_real64 - (dist/real(radius, real64))**2 * 0.7_real64) img(x, y, 1) = int(min(255.0_real64, max(0.0_real64, brightness))) img(x, y, 2) = int(min(255.0_real64, max(0.0_real64, brightness))) img(x, y, 3) = int(min(240.0_real64, max(0.0_real64, brightness * 0.95_real64))) end if end if end do end do ! PPM temporär schreiben open(unit=10, file=filename_ppm, status='replace', action='write', form='formatted') write(10, '(A)') 'P3' write(10, '(I0, 1X, I0)') width, height write(10, '(A)') '255' do y = 1, height do x = 1, width write(10, '(I0, 1X, I0, 1X, I0)') img(x, y, 1), img(x, y, 2), img(x, y, 3) end do end do close(10) ! PPM -> RAW PNG sys_cmd = "pnmtopng " // trim(filename_ppm) // " > " // trim(filename_raw_png) // " && rm " // trim(filename_ppm) call execute_command_line(sys_cmd) ! --- TEXT ZEILENWEISE FORMATIEREN --- write(line_loc, '(A,A,A,F8.5,A,F8.5,A)') & "STANDORT: ", trim(location_name), " (", latitude, " N, ", longitude, " E)" write(line_zeit, '(A,I2.2,A,I2.2,A,I2.2,A,I2.2,A,I2.2,A,I4)') & "ZEIT: ", hour, ":", minute, ":", second, " CEST, ", day, ".", month, ".", year write(line_phase, '(A,A,A,F5.1,A)') & "MONDPHASE: ", trim(phase_str), " (~", phase_pct, "%)" write(line_entf, '(A,I6,A)') & "ENTF.: ~", int(dist_km), " km" write(line_az, '(A,F5.1,A)') & "AZIMUT: ", az, " deg" write(line_alt, '(A,F5.1,A)') & "HOEHE: ", alt, " deg" ! Zeilen zusammenfügen overlay_text = "SIMULATION: MOND-ANSICHT (" // trim(location_name) // ")" // char(10) // char(10) // & trim(line_loc) // char(10) // & trim(line_zeit) // char(10) // & trim(line_phase) // char(10) // & trim(line_entf) // char(10) // & trim(line_az) // char(10) // & trim(line_alt) // char(10) // & "VERGR.: 8x24" // char(10) // "FIELD: 6.5 deg" ! ImageMagick Beschriftung sys_cmd = "convert " // trim(filename_raw_png) // & " -font /usr/share/fonts/truetype/dejavu/DejaVuSansMono.ttf -pointsize 16 -fill white -annotate +30+40 '" // trim(overlay_text) // "' " // & trim(filename_final_png) // " && rm " // trim(filename_raw_png) call execute_command_line(sys_cmd) print *, "Erzeugt: ", trim(filename_final_png) end subroutine render_moon_image ! --- PRÄZISE ASTRONOMISCHE BERECHNUNG (Paul Schlyter / Meeus) --- function get_julian_day_utc(y, m, d, h, min_t, sec_t, utc_offset) result(jd) integer, intent(in) :: y, m, d, h, min_t, sec_t real(real64), intent(in) :: utc_offset real(real64) :: jd, ut integer :: yr, mo, a, b yr = y mo = m if (mo <= 2) then yr = yr - 1 mo = mo + 12 end if a = yr / 100 b = 2 - a + (a / 4) ! Umrechnung der Ortszeit in UTC ut = real(h, real64) + real(min_t, real64)/60.0_real64 + real(sec_t, real64)/3600.0_real64 - utc_offset jd = floor(365.25_real64 * real(yr + 4716, real64)) + & floor(30.6001_real64 * real(mo + 1, real64)) + & real(d, real64) + ut/24.0_real64 + real(b, real64) - 1524.5_real64 end function get_julian_day_utc subroutine calculate_moon_position(jd, lat, lon, alt, az, phase, dist_km, phase_name) real(real64), intent(in) :: jd, lat, lon real(real64), intent(out) :: alt, az, phase, dist_km character(len=*), intent(out) :: phase_name real(real64) :: d, N, inc, w, a, e, M, Ms, ws, Ls, Lm, D_m, F real(real64) :: M_rad, Ms_rad, D_rad, F_rad real(real64) :: d_lambda, lambda, beta, r_earth_radii, obl, obl_rad real(real64) :: lambda_rad, beta_rad, xe, ye, ze, ra, dec real(real64) :: gmst, lmst, ha, lat_rad, dec_rad, ha_rad, d_deg d = jd - 2451543.5_real64 ! Tage seit J2000.0 ! Sonnen-Elemente ws = 282.9404_real64 + 4.70935e-5_real64 * d Ms = mod(356.0470_real64 + 0.9856002585_real64 * d, 360.0_real64) if (Ms < 0.0_real64) Ms = Ms + 360.0_real64 Ls = ws + Ms ! Mond-Bahnelemente N = mod(125.1228_real64 - 0.0529538083_real64 * d, 360.0_real64) inc = 5.1454_real64 w = mod(318.0634_real64 + 0.1643573223_real64 * d, 360.0_real64) a = 60.2665_real64 e = 0.054900_real64 M = mod(115.3654_real64 + 13.0649929509_real64 * d, 360.0_real64) if (M < 0.0_real64) M = M + 360.0_real64 ! Fundamentale Argumente Lm = N + w + M D_m = Lm - Ls F = Lm - N M_rad = M * pi / 180.0_real64 Ms_rad = Ms * pi / 180.0_real64 D_rad = D_m * pi / 180.0_real64 F_rad = F * pi / 180.0_real64 ! Hauptstörungen der ekliptischen Länge (in Grad) d_lambda = -1.274_real64 * sin(M_rad - 2.0_real64*D_rad) + & 0.658_real64 * sin(2.0_real64*D_rad) - & 0.186_real64 * sin(Ms_rad) - & 0.059_real64 * sin(2.0_real64*M_rad - 2.0_real64*D_rad) - & 0.057_real64 * sin(M_rad - 2.0_real64*D_rad + Ms_rad) + & 0.053_real64 * sin(M_rad + 2.0_real64*D_rad) + & 0.046_real64 * sin(2.0_real64*D_rad - Ms_rad) + & 0.041_real64 * sin(M_rad - Ms_rad) - & 0.035_real64 * sin(D_rad) - & 0.031_real64 * sin(M_rad + Ms_rad) lambda = Lm + d_lambda ! Ekliptische Breite beta = 5.128_real64 * sin(F_rad) + 0.280_real64 * sin(M_rad + F_rad) + & 0.277_real64 * sin(M_rad - F_rad) + 0.173_real64 * sin(2.0_real64*D_rad - F_rad) ! Entfernung (in Erdradien, 1 ER = 6378.14 km) r_earth_radii = 60.2665_real64 - 3.300_real64 * cos(M_rad) - & 0.663_real64 * cos(M_rad - 2.0_real64*D_rad) - & 0.370_real64 * cos(2.0_real64*D_rad) dist_km = r_earth_radii * 6378.14_real64 ! Beleuchtung & Phase phase = 50.0_real64 * (1.0_real64 - cos(D_rad)) d_deg = mod(D_m, 360.0_real64) if (d_deg < 0.0_real64) d_deg = d_deg + 360.0_real64 if (phase > 97.0_real64) then phase_name = "VOLLMOND" else if (phase < 3.0_real64) then phase_name = "NEUMOND" else if (d_deg < 180.0_real64) then phase_name = "ZUNEHMEND" else phase_name = "ABNEHMEND" end if ! Ekliptik-Transformation in Äquatorialkoordinaten (RA, DEC) obl = 23.4393_real64 - 0.0000004_real64 * d obl_rad = obl * pi / 180.0_real64 lambda_rad = lambda * pi / 180.0_real64 beta_rad = beta * pi / 180.0_real64 xe = cos(lambda_rad) * cos(beta_rad) ye = sin(lambda_rad) * cos(beta_rad) * cos(obl_rad) - sin(beta_rad) * sin(obl_rad) ze = sin(lambda_rad) * cos(beta_rad) * sin(obl_rad) + sin(beta_rad) * cos(obl_rad) ra = atan2(ye, xe) * 180.0_real64 / pi if (ra < 0.0_real64) ra = ra + 360.0_real64 dec = asin(ze) * 180.0_real64 / pi ! Lokaler Stundenwinkel (GMST -> LMST) gmst = mod(280.46061837_real64 + 360.98564736629_real64 * d, 360.0_real64) lmst = mod(gmst + lon, 360.0_real64) ha = lmst - ra lat_rad = lat * pi / 180.0_real64 dec_rad = dec * pi / 180.0_real64 ha_rad = ha * pi / 180.0_real64 ! Höhe (Altitude) & Azimut alt = asin(sin(lat_rad)*sin(dec_rad) + cos(lat_rad)*cos(dec_rad)*cos(ha_rad)) * 180.0_real64 / pi az = atan2(-sin(ha_rad), tan(dec_rad)*cos(lat_rad) - sin(lat_rad)*cos(ha_rad)) * 180.0_real64 / pi if (az < 0.0_real64) az = az + 360.0_real64 end subroutine calculate_moon_position pure function lowercase(str) result(res) character(len=*), intent(in) :: str character(len=len(str)) :: res integer :: j, code do j = 1, len(str) code = iachar(str(j:j)) if (code >= 65 .and. code <= 90) then res(j:j) = achar(code + 32) else res(j:j) = str(j:j) end if end do end function lowercase end program generate_moon