doomsday.f
1 ! @section Synopsis
2 !
3 ! Doomsday algorithm giving the day of the week for a date.
4 !
5 ! @author J. Marcel van der Veer
6 !
7 ! @section copyright
8 !
9 ! This file is part of VIF - vintage fortran compiler.
10 ! Copyright 2020-2025 J. Marcel van der Veer <algol68g@xs4all.nl>.
11 !
12 ! @section license
13 !
14 ! This program is free software; you can redistribute it and/or modify it
15 ! under the terms of the gnu general public license as published by the
16 ! free software foundation; either version 3 of the license, or
17 ! (at your option) any later version.
18 !
19 ! This program is distributed in the hope that it will be useful, but
20 ! without any warranty; without even the implied warranty of merchantability
21 ! or fitness for a particular purpose. See the GNU general public license for
22 ! more details. You should have received a copy of the GNU general public
23 ! license along with this program. If not, see <http://www.gnu.org/licenses/>.
24 !
25 subroutine doomsday_gregorian ( y, m, d, w )
26
27 c DOOMSDAY_GREGORIAN(): weekday given any date in Gregorian calendar.
28 c
29 c Discussion:
30 c
31 c This procedure does not include any procedure to switch to the Julian
32 c calendar for dates early enough that that calendar was used instead.
33 c
34 c Licensing:
35 c
36 c This code is distributed under the MIT license.
37 c
38 c Modified:
39 c
40 c 28 May 2012
41 c
42 c Author:
43 c
44 c John Burkardt
45 c
46 c Reference:
47 c
48 c John Conway,
49 c Tomorrow is the Day After Doomsday,
50 c Eureka,
51 c Volume 36, October 1973, pages 28-31.
52 c
53 c Parameters:
54 c
55 c Input, integer Y, M, D, the year, month and day of the date.
56 c Note that the year must be positive.
57 c
58 c Output, integer W, the weekday of the date.
59 c
60 implicit none
61
62 integer anchor(4)
63 integer c
64 integer d
65 integer drd
66 integer drdr
67 integer i4_wrap
68 integer l
69 integer m
70 integer mdoom(12)
71 integer w
72 integer y
73 integer ydoom
74 logical year_is_leap_gregorian
75 integer yy
76 integer yy12d
77 integer yy12r
78 integer yy12r4d
79
80 save anchor
81 save mdoom
82
83 data anchor / 1, 6, 4, 3 /
84 data mdoom / 3, 28, 0, 4, 9, 6, 11, 8, 5, 10, 7, 12 /
85 c
86 c Refuse to handle Y <= 0.
87 c
88 if ( y .le. 0 ) then
89 write ( *, '(a)' ) ' '
90 write ( *, '(a)' ) 'DOOMSDAY_GREGORIAN - Fatal error!'
91 write ( *, '(a)' ) ' Y <= 0.'
92 stop
93 end if
94 c
95 c Determine the century C.
96 c
97 c = y / 100
98 c
99 c Determine the last two digits of the year, YY
100 c
101 yy = mod ( y, 100 )
102 c
103 c Divide the last two digits of the year by 12.
104 c
105 yy12d = yy / 12
106 yy12r = mod ( yy, 12 )
107 yy12r4d = yy12r / 4
108 drd = yy12d + yy12r + yy12r4d
109
110 drdr = mod ( drd, 7 )
111 ydoom = anchor( mod ( c-1, 4 ) + 1 ) + drdr
112 ydoom = i4_wrap ( ydoom, 1, 7 )
113 c
114 c If M = 1 or 2, and Y is a leap year, add 1.
115 c
116 if ( ( m .eq. 1 .or. m .eq. 2 ) .and.
117 & year_is_leap_gregorian ( y ) ) then
118 l = 1
119 else
120 l = 0
121 end if
122
123 w = ydoom + ( d - mdoom(m) - l )
124 w = i4_wrap ( w, 1, 7 )
125
126
127 return
128 end
129
130 function i4_modp ( i, j )
131
132 cc I4_MODP returns the nonnegative remainder of integer division.
133 c
134 c Discussion:
135 c
136 c If
137 c NREM = I4_MODP ( I, J )
138 c NMULT = ( I - NREM ) / J
139 c then
140 c I = J * NMULT + NREM
141 c where NREM is always nonnegative.
142 c
143 c The MOD function computes a result with the same sign as the
144 c quantity being divided. Thus, suppose you had an angle A,
145 c and you wanted to ensure that it was between 0 and 360.
146 c Then mod(A,360) would do, if A was positive, but if A
147 c was negative, your result would be between -360 and 0.
148 c
149 c On the other hand, I4_MODP(A,360) is between 0 and 360, always.
150 c
151 c Example:
152 c
153 c I J MOD I4_MODP Factorization
154 c
155 c 107 50 7 7 107 = 2 * 50 + 7
156 c 107 -50 7 7 107 = -2 * -50 + 7
157 c -107 50 -7 43 -107 = -3 * 50 + 43
158 c -107 -50 -7 43 -107 = 3 * -50 + 43
159 c
160 c Licensing:
161 c
162 c This code is distributed under the MIT license.
163 c
164 c Modified:
165 c
166 c 30 December 2006
167 c
168 c Author:
169 c
170 c John Burkardt
171 c
172 c Parameters:
173 c
174 c Input, integer I, the number to be divided.
175 c
176 c Input, integer J, the number that divides I.
177 c
178 c Output, integer I4_MODP, the nonnegative remainder when I is
179 c divided by J.
180 c
181 implicit none
182
183 integer i
184 integer i4_modp
185 integer j
186 integer value
187
188 if ( j .eq. 0 ) then
189 call xerabt ('i4_modp: illegal divisor', 1)
190 end if
191
192 value = mod ( i, j )
193
194 if ( value .lt. 0 ) then
195 value = value + abs ( j )
196 end if
197
198 i4_modp = value
199
200 return
201 end
202
203 function i4_wrap ( ival, ilo, ihi )
204
205 c I4_WRAP forces an I4 to lie between given limits by wrapping.
206 c
207 c Example:
208 c
209 c ILO = 4, IHI = 8
210 c
211 c I Value
212 c
213 c -2 8
214 c -1 4
215 c 0 5
216 c 1 6
217 c 2 7
218 c 3 8
219 c 4 4
220 c 5 5
221 c 6 6
222 c 7 7
223 c 8 8
224 c 9 4
225 c 10 5
226 c 11 6
227 c 12 7
228 c 13 8
229 c 14 4
230 c
231 c Licensing:
232 c
233 c This code is distributed under the MIT license.
234 c
235 c Modified:
236 c
237 c 31 December 2006
238 c
239 c Author:
240 c
241 c John Burkardt
242 c
243 c Parameters:
244 c
245 c Input, integer IVAL, an integer value.
246 c
247 c Input, integer ILO, IHI, the desired bounds for the integer value.
248 c
249 c Output, integer I4_WRAP, a "wrapped" version of IVAL.
250 c
251 implicit none
252
253 integer i4_modp
254 integer i4_wrap
255 integer ihi
256 integer ilo
257 integer ival
258 integer jhi
259 integer jlo
260 integer value
261 integer wide
262
263 jlo = min ( ilo, ihi )
264 jhi = max ( ilo, ihi )
265
266 wide = jhi - jlo + 1
267
268 if ( wide .eq. 1 ) then
269 value = jlo
270 else
271 value = jlo + i4_modp ( ival - jlo, wide )
272 end if
273
274 i4_wrap = value
275
276 return
277 end
278
279 subroutine weekday_to_name_common ( w, s )
280
281 c WEEKDAY_TO_NAME_COMMON returns the name of a Common weekday.
282 c
283 c Licensing:
284 c
285 c This code is distributed under the MIT license.
286 c
287 c Modified:
288 c
289 c 28 August 1999
290 c
291 c Author:
292 c
293 c John Burkardt
294 c
295 c Parameters:
296 c
297 c Input, integer W, the weekday index.
298 c
299 c Output, character * ( * ) S, the weekday name.
300 c
301 implicit none
302
303 character * 9 name(7)
304 character * ( * ) s
305 integer w
306 integer w2
307
308 save name
309
310 data name /
311 & 'Sunday ', 'Monday ', 'Tuesday ', 'Wednesday',
312 & 'Thursday ', 'Friday ', 'Saturday ' /
313 c
314 c Make a local copy of the weekday number.
315 c
316 w2 = w
317 c
318 c Return the weekday name.
319 c
320 s = name ( w2 )
321
322 return
323 end
324
325 subroutine weekday_values ( n_data, y, m, d, w )
326
327 c WEEKDAY_VALUES returns the day of the week for various dates.
328 c
329 c Discussion:
330 c
331 c The CE or Common Era calendar is used, under the
332 c hybrid Julian/Gregorian Calendar, with a transition from Julian
333 c to Gregorian. The day after 04 October 1582 was 15 October 1582.
334 c
335 c The year before 1 AD or CE is 1 BC or BCE. In this data set,
336 c years BC/BCE are indicated by a negative year value.
337 c
338 c Licensing:
339 c
340 c This code is distributed under the MIT license.
341 c
342 c Modified:
343 c
344 c 26 May 2012
345 c
346 c Author:
347 c
348 c John Burkardt
349 c
350 c Reference:
351 c
352 c Edward Reingold, Nachum Dershowitz,
353 c Calendrical Calculations: The Millennium Edition,
354 c Cambridge University Press, 2001,
355 c ISBN: 0 521 77752 6
356 c LC: CE12.R45.
357 c
358 c Parameters:
359 c
360 c Input/output, integer N_DATA. The user sets N_DATA to 0
361 c before the first call. On each call, the routine increments N_DATA by 1,
362 c and returns the corresponding data; when there is no more data, the
363 c output value of N_DATA will be 0 again.
364 c
365 c Output, integer Y, M, D, the Common Era date.
366 c
367 c Output, integer W, the day of the week. Sunday = 1.
368 c
369 implicit none
370
371 integer n_max
372 parameter ( n_max = 34 )
373
374 integer d
375 integer d_vec(n_max)
376 integer m
377 integer m_vec(n_max)
378 integer n_data
379 integer w
380 integer w_vec(n_max)
381 integer y
382 integer y_vec(n_max)
383
384 save d_vec
385 save m_vec
386 save w_vec
387 save y_vec
388
389 data d_vec /
390 & 30,
391 & 8,
392 & 26,
393 & 3,
394 & 7,
395 & 18,
396 & 7,
397 & 19,
398 & 14,
399 & 18,
400 & 16,
401 & 3,
402 & 26,
403 & 20,
404 & 4,
405 & 25,
406 & 31,
407 & 9,
408 & 24,
409 & 10,
410 & 30,
411 & 24,
412 & 19,
413 & 2,
414 & 27,
415 & 19,
416 & 25,
417 & 29,
418 & 19,
419 & 7,
420 & 17,
421 & 25,
422 & 10,
423 & 18 /
424 data m_vec /
425 & 7,
426 & 12,
427 & 9,
428 & 10,
429 & 1,
430 & 5,
431 & 11,
432 & 4,
433 & 10,
434 & 5,
435 & 3,
436 & 3,
437 & 3,
438 & 4,
439 & 6,
440 & 1,
441 & 3,
442 & 9,
443 & 2,
444 & 6,
445 & 6,
446 & 7,
447 & 6,
448 & 8,
449 & 3,
450 & 4,
451 & 8,
452 & 9,
453 & 4,
454 & 10,
455 & 3,
456 & 2,
457 & 11,
458 & 7 /
459 data w_vec /
460 & 1,
461 & 4,
462 & 4,
463 & 1,
464 & 4,
465 & 2,
466 & 7,
467 & 1,
468 & 7,
469 & 1,
470 & 6,
471 & 7,
472 & 6,
473 & 1,
474 & 1,
475 & 4,
476 & 7,
477 & 7,
478 & 7,
479 & 4,
480 & 1,
481 & 6,
482 & 1,
483 & 2,
484 & 4,
485 & 1,
486 & 1,
487 & 2,
488 & 2,
489 & 5,
490 & 3,
491 & 1,
492 & 4,
493 & 1 /
494 data y_vec /
495 & - 587,
496 & - 169,
497 & 70,
498 & 135,
499 & 470,
500 & 576,
501 & 694,
502 & 1013,
503 & 1066,
504 & 1096,
505 & 1190,
506 & 1240,
507 & 1288,
508 & 1298,
509 & 1391,
510 & 1436,
511 & 1492,
512 & 1553,
513 & 1560,
514 & 1648,
515 & 1680,
516 & 1716,
517 & 1768,
518 & 1819,
519 & 1839,
520 & 1903,
521 & 1929,
522 & 1941,
523 & 1943,
524 & 1943,
525 & 1992,
526 & 1996,
527 & 2038,
528 & 2094 /
529
530 if ( n_data .lt. 0 ) then
531 n_data = 0
532 end if
533
534 n_data = n_data + 1
535
536 if ( n_max .lt. n_data ) then
537 n_data = 0
538 y = 0
539 m = 0
540 d = 0
541 w = 0
542 else
543 y = y_vec(n_data)
544 m = m_vec(n_data)
545 d = d_vec(n_data)
546 w = w_vec(n_data)
547 end if
548
549 return
550 end
551
552 function year_is_leap_gregorian ( y )
553
554 c YEAR_IS_LEAP_GREGORIAN returns TRUE if the Gregorian year was a leap year.
555 c
556 c Licensing:
557 c
558 c This code is distributed under the MIT license.
559 c
560 c Modified:
561 c
562 c 21 May 2000
563 c
564 c Author:
565 c
566 c John Burkardt
567 c
568 c Parameters:
569 c
570 c Input, integer Y, the year to be checked.
571 c
572 c Output, logical YEAR_IS_LEAP_GREGORIAN, TRUE if the year was a leap year,
573 c FALSE otherwise.
574 c
575 implicit none
576
577 integer y
578 logical year_is_leap_gregorian
579
580 if ( y .le. 0 ) then
581 call xerabt ('year_is_leap_gregorian: invalid year', 1)
582 end if
583
584 if ( mod ( y, 400 ) .eq. 0 ) then
585 year_is_leap_gregorian = .true.
586 else if ( mod ( y, 100 ) .eq. 0 ) then
587 year_is_leap_gregorian = .false.
588 else if ( mod ( y, 4 ) .eq. 0 ) then
589 year_is_leap_gregorian = .true.
590 else
591 year_is_leap_gregorian = .false.
592 end if
593
594 return
595 end
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
