年月日与年积日（日序数）互转——今夕何夕（今天是本年度第几天）？

hillside

    “今夕何夕？”、“今天是本年度第几天？”似乎是追星族们推算命运星座等的游戏之谈。其实，它在GPS领域以及气候、水文等相关地学领域已有重要的应用。
     现附上一个国外程序（国内用户改动了输出显示）以供有意者阅读与使用，摘录如下：

      program doy
*     This is generic program which will read a varierty of
*     date types from the runstring and return the other date
*     types.

*   date(5)     - Calender date
*   day_of_year - Day of year
*   rcpar       - Get runstring parameters
*   len_run     - Length of string entry
*   ierr        - IOSTAT error flag
*   num_run     - Number of arguments in the runstring
*   i           - Loop counter
*   gps_week    - GPS week number starts with week zero on 1980/1/5 Midnight
*   gps_dow     - GPS doy of week
*   gps_sow     - GPS seconds into week
*   gps_start_mjd - MJD of the start of the GPS week count (80/1/5)
*   indx        - Postition of character in string
c   Add Decimal Year in output, K. Feigl May 97
c   MOD TAH 970520: Changed gps_dow to run from 0-6 (instead of 1-7).
c   MOD SCM 971029: Changed some of the math from REAL*8 to INTEGER*4 for LINUX.
c   MOD SCM 971115: Fixed more problems with mixing REAL and INTEGER arithmatic.

      integer*4 date(5), day_of_year, rcpar, len_run, ierr, num_run, i,
     .          gps_week, gps_dow, gps_sow, gps_start_mjd , indx

*   jd          - Julian date
*   mjd         - Modified Julian date.
*   sectag      - Seconds tag

      real*8 jd, mjd, sectag,dyear

*   runstring(5)    - Five elments of the string.

      character*40 end_88
      character*40 runstring(5)
      character*4 day_of_week(7)
      data gps_start_mjd / 44243 /
      data day_of_week / 'Sun ', 'Mon ', 'Tue ', 'Wed ',
     .                   'Thu ', 'Fri ', 'Sat '  /
      write(*,*) '本功能实现 年月日<－>年积日 －> GPS周和小数点年 '
      write(*,*) '请输入 相应日期 <－>  '
write(*,*) '1＝“年 月 日”'
write(*,*) '2＝“年 年积日”'
......

气象家园相关帖子：
[经验总结] 日序数的求算公式

[源代码] ~闲来无聊，写个求每个月第一天是一年中第几天的小程序

附：一篇介绍在GPS时（GPS周和周内秒）和公历日期之间进行转换的文章（注：下文中“我”指写作下文的原作者。文中所提多个MATLAB程序见于附件）
                                                                                    “1 概述
        最近在学习GPS解算算法时需要在GPS时（GPS周和周内秒）和公历日期之间进行转换，于是就整理了一些时间转换的小程序。本文介绍了GPS时、公历、儒略日（JD）、简化儒略日（MJD）之间的转换函数。
gps2cal 将GPS周和周内秒转换到公历时间
gps2cal1 由公历日期和GPS周内秒计算公历时间
cal2gps 将公历GPS时间转换到GPS周和周内秒
mjd2cal 将简化儒略日转换到公历时间
cal2mjd 将公历时间转换到简化儒略日
jd2cal 从儒略日计算公历时间
cal2jd 将公历时间转换到儒略日
cal2wd 借助MJD，由公历年月日推算星期几，按照格里高利十三世的历法改革去掉1582年10月5日至14日
cal2wd1 由公历年月日推算星期几，按照英国人的做法去掉1752年9月3日至13日
cal2wd2 不借助MJD，由公历年月日推算星期几，去掉1582年10月5日至14日

        在MJD基础上推算星期是很简单的，我顺带写出了从公历推算星期的函数cal2wd。不借助MJD也可以推算星期，例如安竹林同学在《再谈星期的计算》（《程序员》2001年第4期）介绍的方法。本文给出了这个方法的matlab版本cal2wd1。 cal2wd1按照英国人的历法去掉了公历中的公元1752年9月3日至13日。我又写了cal2wd2函数，这个函数按照格里高利十三世的做法去掉1582年10月5日至14日。
        我在2005年写过一篇《时标与历法》，介绍过一些历法的基础知识。本文假设读者了解该文介绍的一些概念，例如原子时标、UTC、闰秒、儒略历、儒略日、简化儒略日等。

                                                                                                   2 GPS周和周内秒

        GPS时从1980年1月6日0时开始计时。GPS时的秒长度与UTC一样，采用原子时标。但GPS时是连续的，不调整闰秒。在GPS数据中，GPS时通常被表示为GPS周和周内秒。例如：
>> gpst=cal2gps([2009 3 13 22 38 10])
gpst =
        1522      513490
>> cal=gps2cal(gpst)
cal =
        2009           3          13          22          38          10
        即2009年3月13日22:38:10对应gps周1522，周内秒513490。 GPS数据为什么要采用周和周内秒的表示方式呢？
        因为这样便于存储和处理。GPS系统内部每隔1.5s会产生一个叫作X1历元的定时事件。系统通过对X1历元计数的方式计时，这种计时方式被称作Z计数。 Z计数有29bit，高10bit是周数，低19bit是周内秒。一个星期有604800秒，19bit所能表示的最长时间是(2^19-1)*1.5=786430.5秒，够用啦。 GPS卫星和接收机之间通过Z计数传递时间，所以GPS数据文件中会用到GPS周和周内秒。
我在做数据解算时，有时知道观测日期和周内秒，所以写了一个由公历日期和GPS周内秒计算公历时间的函数gps2cal1：
>> cal=gps2cal1([2009 3 13],513490)
cal =
        2009           3          13          22          38          10
        请注意这3个函数里的公历时间没有调整闰秒。在GPS数据（例如RINEX格式）中，公历时间也不调整闰秒。

                                                                                                 3 儒略日和简化儒略日
        我在《时标与历法》中介绍过儒略历、儒略日（JD）和简化儒略日（MJD），这里就不再重复。只是简单说明一下公历的不连续性。公元1582年3月1日，罗马教皇格里高利十三世在颁布格里高利历（即我们现在使用的公历）时，为了消除已经存在的误差，宣布去掉1582年10月5日至14日。即1582年10月4日的23:59的下一秒是1582年10月15日的00:00。
       我们通常对1582年10月15日0时之后的时间使用格里高利历，对1582年10月4日的23:59之前的时间使用儒略历。这里的儒略历是指屋大维改革过的儒略历，即奥古斯都历。它和格里高利历除了置闰方式不同，其它是完全一样的。儒略历是4年一闰，400年100闰。公历是400年97闰，扣掉了3个不能被400整除的世纪年。
        例如：香港天文台说儒略日“以倒推到公元前4713年1月1日格林尼治平时正午为起算日期”。这里的公元前4713年1月1日就是按照儒略历倒推的日期。但没有必要特殊说明儒略历。本文提供的mjd2cal、cal2mjd、jd2cal、cal2jd函数都使用了规则：公元1582年10月4日24:00点之前使用儒略历，公元1582年10月15日00:00点之后使用格里高利历。下面是一些使用示例：
>> jd2cal(0)
ans =
       -4713           1           1          12           0           0
>> mjd2cal(0)
ans =
        1858          11          17           0           0           0
>> cal2jd([1582 10 4 12 0 0])
ans =
     2299160
>> cal2jd([1582 10 15 12 0 0])
ans =
     2299161
>> cal2mjd([1582 10 4 0 0 0])
ans =
     -100841
>> cal2mjd([1582 10 15 0 0 0])
ans =
     -100840
        网上有很多介绍闰年的文章会以公元1000年为例，说它不是闰年。事实上，按照1582年10月4日24:00点之前使用儒略历的规则，公元1000年是闰年。也就是说公元1000年2月29日在历史上是存在的。例如：
>> cal2jd([1000 2 28 12 0 0])
ans =
     2086366
>> cal2jd([1000 2 29 12 0 0])
ans =
     2086367
>> cal2jd([1000 3 1 12 0 0])
ans =
     2086368
3 星期推算
在MJD的基础上，我们可以很方便地推算星期。函数cal2wd只需要两行代码：
mjd=floor(cal2mjd(cal));
% 2009年3月9日（MJD 54899）是星期一
wd=mod(mjd-54899,7)+1;
因为在JD或MJD中，1582年10月5日至14日是不存在的。所以1582年10月15日是星期五，1582年10月4日是星期四：
>> cal2wd([1582 10 15])
ans =
     5
>> cal2wd([1582 10 4])
ans =
     4
        不借助MJD，也可以直接推算星期。安竹林同学在《再谈星期的计算》一文中介绍了一个简单易行的方法。我将其写成matlab函数cal2wd1。不过在安同学的代码中，1582年10月5日至14日是存在的，1752年9月3日至13日是不存在的，这是怎么回事呢？
        事情是这样的。1582年3月1日格里高利十三世颁布格里高利历后，英国人认为从历史上去掉10天是很荒唐的事情，没有接受。其实更合理的解释是：格里高利历是罗马教皇颁布的，信奉天主教的意大利、波兰、西班牙、葡萄牙都国家都很痛快地接受了。但英国是新教国家，所以不愿意接受。直到1752年，英国人才决定采用格里高利历，因为这时儒略历又多算了一天，所以英国人从历史上去掉了11天。在英国人看来，1752年9月2日的次日是1752年9月14日。
        在我看来，还是去掉1582年10月5日至14日更常用一些。所以又在cal2wd1基础上写了cal2wd2。 cal2wd2和cal2wd功能相同，区别仅是前者不依赖cal2mjd函数。cal2wd和cal2wd1也就是对1582年10月5日至1752年9月13日的推算结果不同。它们对其它日期的推算结果都是相同的，例如它们推算的公元1年1月1日都是星期六：
>> cal2wd([1752 9 14])
ans =
     4
>> cal2wd1([1752 9 14])
ans =
     4
>> cal2wd([1 1 1])
ans =
     6
>> cal2wd1([1 1 1])
ans =
     6
                                                                                       4 结束语
        本文提到的源代码可以从这里下载。其中最复杂的cal2jd和jd2cal是根据www.moshier.net上的代码修改得到。其实我起初只是想写两个转换gps时间的函数。好奇心让我在半夜查阅相关资料，并用本文记录了这些函数的一些背景。
本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/fmddlmyy/archive/2009/03/15/3991315.aspx”

说明：本帖后面的内容可能更适合放在3S版，但因年积日转换程序为FORTRAN版，故作为内容交叉帖置放于此。

                                                                                     doy.exe用法介绍



                                                                   年积日查询表
                                        http://amsu.cira.colostate.edu/julian.html              
         以下为非闰年摘录，欲知年积日对照表全部，请点击上述网址。

liutaoz21

感谢楼主分享！

infocert

感谢楼主分享！

kongfeng0824

很详细呀！！！！！！！！！！！！！！！！

mofangbao

很不错的功能，这个很容易会用到

cuhum

有意思！

shine090

好可爱又好复杂

醉卧不知日已暮

很不错，用得着

weiwei110

牛人啊  

zane

  赞！