地球仪8字曲线估算时间方程式公式差

!***************************************************************
!***************************************************************
! 地球仪8字曲线估算时间方程式公式差
!***************************************************************
! SINOMACH
!
!地球仪8字曲线估算时间方程式公式差于一致的24小时和太阳的实际位置,
!形成数据用于图形, 点数可上至五十万五千（505，000）.
!
! 周勇 20221125
!***************************************************************
!c $$$ Run Output:
!***************************************************************
!***************************************************************
Program main

!*********************************************************************72
!
!c MAIN is the main program for ANALEMMA.
!
! Local parameters:
!
! Local, double precision ECC, the orbital eccentricity.
!
! Local, double precision LON, the longitude of the perihelion in radians.
!
! Local, double precision OBLIQ, the obliquity in radians.
!
Implicit None

Character *(255) arg
Character *(80) command_filename
Integer command_unit
Character *(80) data_filename
Integer data_unit
Double Precision days
Parameter (days=365.242D+00)
Double Precision dec
Double Precision degrees
Parameter (degrees=3.141592653589793D+00/180.0D+00)
Double Precision ecc
Double Precision eot
Double Precision f
Integer i
Double Precision lon
Integer n
Integer numarg
Double Precision obliq
Double Precision pi
Parameter (pi=3.141592653589793D+00)
Double Precision t
Double Precision tau
Double Precision theta
Double Precision x1
Double Precision x2
Double Precision x3
Double Precision y1
Double Precision y2
Double Precision y3
Double Precision z1
Double Precision z2
Double Precision z3

Call timestamp()
Write (, '(a)') ''
Write (, '(a)') 'ANALEMMA'
Write (, '(a)') ' FORTRAN95, 2013&2018 version'
Write (, '(a)') ' Compute and plot the analemma, equation of time,'
Write (, '(a)') ' and declination.'
Write (, '(a)') ' This program is d on a C program by Brian Tung.'
!
! Parse the arguments
!
numarg = iargc()

If (numarg<1) Then
ecc = 0.01671D+00
Else
i = 1
Call getarg(i, arg)
Call s_to_r8(arg, ecc)
End If

If (numarg<2) Then
lon = 1.347D+00
Else
i = 2
Call getarg(i, arg)
Call s_to_r8(arg, lon)
lon = lon*pi/180.0D+00
End If

If (numarg<3) Then
obliq = 0.4091D+00
Else
i = 3
Call getarg(i, arg)
Call s_to_r8(arg, obliq)
obliq = obliq*pi/180.0D+00
End If
!
! Compute the data.
!
data_filename = 'analemma_data.txt'
Call get_unit(data_unit)
Open (Unit=data_unit, File=data_filename)

n = 25000

Do i = 0, n - 1

f = dble(i)/dble(n)
tau = 2.0D+00*pi*f

!
! Set theta to the current longitude.
!
theta = atan2(sqrt(1.0D+00-ecc*ecc)*sin(tau), cos(tau)-ecc)
!
! Rotate clockwise in XY plane by theta, corrected by lon.
!
x1 = cos(theta-(lon-pi/2.0D+00))
y1 = sin(theta-(lon-pi/2.0D+00))
z1 = 0.0D+00
!
! Rotate counter-clockwise in XZ plane by obliq.
!
x2 = cos(obliq)*x1 + sin(obliq)*z1
y2 = y1
z2 = -sin(obliq)*x1 + cos(obliq)z1
!
! Set t equal to real time from tau and
! rotate counter-clockwise by t, corrected by lon
!
t = tau - eccsin(tau)
x3 = cos(t-(lon-pi/2.0D+00))*x2 + sin(t-(lon-pi/2.0D+00))*y2
y3 = -sin(t-(lon-pi/2.0D+00))*x2 + cos(t-(lon-pi/2.0D+00))*y2
z3 = z2

eot = -atan2(y3, x3)*4.0D+00/degrees*days/(days+1.0D+00)
dec = asin(z3)/degrees

!
! Print results in minutes early/late and degrees north/south
!
Write (data_unit, '(2x,g15.9,2x,g15.9,2x,g15.9)') t/(2.0D+00*pi), eot, dec

End Do

Close (Unit=data_unit)

Write (, '(a)') ''
Write (, '(a)') ' Created data file "' // trim(data_filename) // '".'
!
! Create the command file.
!
command_filename = 'analemma_commands.txt'
Call get_unit(command_unit)
Open (Unit=command_unit, File=command_filename)
Write (command_unit, '(a)') 'set term png'
Write (command_unit, '(a)') 'set grid'
Write (command_unit, '(a)') 'set style data lines'
Write (command_unit, '(a)') 'unset key'
Write (command_unit, '(a)') 'set output "eot.png"'
Write (command_unit, '(a)') 'set xlabel "<---Normalized Date--->"'
Write (command_unit, '(a)') 'set ylabel "<---Minutes Early/Late--->"'
Write (command_unit, '(a)') 'set "The equation of time"'
Write (command_unit, '(a)') 'plot "' // trim(data_filename) //&
'" using 1:2 with lines'
Write (command_unit, '(a)') 'set output "declination.png"'
Write (command_unit, '(a)') 'set xlabel "<---Normalized Date--->"'
Write (command_unit, '(a)') 'set ylabel "<---Degrees North/South--->"'
Write (command_unit, '(a)') 'set "Declination"'
Write (command_unit, '(a)') 'plot "' // trim(data_filename) // &
'" using 1:3 with lines'
Write (command_unit, '(a)') 'set output "analemma.png"'
Write (command_unit, '(a)') 'set xlabel "<---Minutes Early/Late--->"'
Write (command_unit, '(a)') 'set ylabel "<---Degrees North/South--->"'
Write (command_unit, '(a)') 'set "The analemma"'
Write (command_unit, '(a)') 'plot "' // trim(data_filename) //&
'" using 2:3 with lines'
Write (command_unit, '(a)') 'quit'

Close (Unit=command_unit)

Write (, '(a)') ' Created command file "' // trim(command_filename) // '".'
!
! Terminate.
!
Write (, '(a)') ''
Write (, '(a)') 'ANALEMMA'
Write (, '(a)') ' Normal end of execution.'
Write (*, '(a)') ''
Call timestamp()

Stop
End Program main
Subroutine get_unit(iunit)

!*********************************************************************72
!
!c GET_UNIT returns a free FORTRAN unit number.
!
! Discussion:
!
! A "free" FORTRAN unit number is a value between 1 and 99 which
! is not currently associated with an I/O device. A free FORTRAN unit
! number is needed in order to open a file with the OPEN command.
!
! If IUNIT = 0, then no free FORTRAN unit could be found, although
! all 99 units were checked (except for units 5, 6 and 9, which
! are commonly reserved for console I/O).
!
! Otherwise, IUNIT is a value between 1 and 99, representing a
! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6
! are special, and will never return those values.
!
! Parameters:
!
! Output, integer IUNIT, the free unit number.
!
Implicit None

Integer i
Integer iunit
Logical value

iunit = 0

Do i = 1, 99

If (i/=5 .And. i/=6 .And. i/=9) Then

  Inquire (Unit=i, Opened=value, Err=10)

  If (.Not. value) Then
    iunit = i
    Return
  End If

End If

10 Continue

End Do

Return
End Subroutine get_unit
Subroutine s_to_r8(s, r8)

!*********************************************************************72
!
!c S_TO_R8 reads an R8 value from a string.
!
! Discussion:
!
! An "R8" value is simply a real number to be stored as a
! variable of type "double precision".
!
! The routine will read as many characters as possible until it reaches
! the end of the string, or encounters a character which cannot be
! part of the number.
!
! Legal input is:
!
! 1 blanks,
! 2 '+' or '-' sign,
! 2.5 blanks
! 3 integer part,
! 4 decimal point,
! 5 fraction part,
! 6 'E' or 'e' or 'D' or 'd', exponent marker,
! 7 exponent sign,
! 8 exponent integer part,
! 9 exponent decimal point,
! 10 exponent fraction part,
! 11 blanks,
! 12 final comma or semicolon,
!
! with most quantities optional.
!
! Example:
!
! S R8
!
! '1' 1.0
! ' 1 ' 1.0
! '1A' 1.0
! '12,34,56' 12.0
! ' 34 7' 34.0
! '-1E2ABCD' -100.0
! '-1X2ABCD' -1.0
! ' 2E-1' 0.2
! '23.45' 23.45
! '-4.2E+2' -420.0
! '17d2' 1700.0
! '-14e-2' -0.14
! 'e2' 100.0
! '-12.73e-9.23' -12.73 * 10.0^(-9.23)
!
! Parameters:
!
! Input, character * ( * ) S, the string containing the
! data to be read. Reading will begin at position 1 and
! terminate at the end of the string, or when no more
! characters can be read to form a legal real. Blanks,
! commas, or other nonnumeric data will, in particular,
! cause the conversion to halt.
!
! Output, double precision R8, the value read from the string.
!
Implicit None

Character c
Integer ierror
Integer ihave
Integer isgn
Integer iterm
Integer jbot
Integer jsgn
Integer jtop
Integer length
Integer ndig
Double Precision r8
Double Precision rbot
Double Precision rexp
Double Precision rtop
Character () s
Integer s_length
Character tab
Parameter (tab=char(9))

s_length = len_trim(s)

ierror = 0
r8 = 0.0D+00
length = -1
isgn = 1
rtop = 0
rbot = 1
jsgn = 1
jtop = 0
jbot = 1
ihave = 1
iterm = 0

10 Continue

length = length + 1

If (s_length<length+1) Then
Goto 20
End If

c = s(length+1:length+1)
!
! Blank character.
!
If (c==' ' .Or. c==tab) Then

If (ihave==2) Then

Else If (ihave==6 .Or. ihave==7) Then
  iterm = 1
Else If (1<ihave) Then
  ihave = 11
End If

!
! Comma.
!
Else If (c==',' .Or. c==';') Then

If (ihave/=1) Then
  iterm = 1
  ihave = 12
  length = length + 1
End If

!
! Minus sign.
!
Else If (c=='-') Then

If (ihave==1) Then
  ihave = 2
  isgn = -1
Else If (ihave==6) Then
  ihave = 7
  jsgn = -1
Else
  iterm = 1
End If

!
! Plus sign.
!
Else If (c=='+') Then

If (ihave==1) Then
  ihave = 2
Else If (ihave==6) Then
  ihave = 7
Else
  iterm = 1
End If

!
! Decimal point.
!
Else If (c=='.') Then

If (ihave<4) Then
  ihave = 4
Else If (6<=ihave .And. ihave<=8) Then
  ihave = 9
Else
  iterm = 1
End If

!
! Scientific notation exponent marker.
!
Else If (c=='E' .Or. c=='e' .Or. c=='D' .Or. c=='d') Then

If (ihave<6) Then
  ihave = 6
Else
  iterm = 1
End If

!
! Digit.
!
Else If (ihave<11 .And. lle('0',c) .And. lle(c,'9')) Then

If (ihave<=2) Then
  ihave = 3
Else If (ihave==4) Then
  ihave = 5
Else If (ihave==6 .Or. ihave==7) Then
  ihave = 8
Else If (ihave==9) Then
  ihave = 10
End If

ndig = ichar(c) - 48

If (ihave==3) Then
  rtop = 10.0D+00*rtop + dble(ndig)
Else If (ihave==5) Then
  rtop = 10.0D+00*rtop + dble(ndig)
  rbot = 10.0D+00*rbot
Else If (ihave==8) Then
  jtop = 10*jtop + ndig
Else If (ihave==10) Then
  jtop = 10*jtop + ndig
  jbot = 10*jbot
End If

!
! Anything else is regarded as a terminator.
!
Else
iterm = 1
End If
!
! If we haven't seen a terminator, and we haven't examined the
! entire string, go get the next character.
!
If (iterm==1) Then
Goto 20
End If

Goto 10

20 Continue
!
! If we haven't seen a terminator, and we have examined the
! entire string, then we're done, and LENGTH is equal to S_LENGTH.
!
If (iterm/=1 .And. length+1==s_length) Then
length = s_length
End If
!
! Number seems to have terminated. Have we got a legal number?
! Not if we terminated in states 1, 2, 6 or 7!
!
If (ihave==1 .Or. ihave==2 .Or. ihave==6 .Or. ihave==7) Then
ierror = ihave
Write (*, '(a)') ' '
Write (*, '(a)') 'S_TO_R8 - Serious errorc'
Write (*, '(a)') ' Illegal or nonnumeric input:'
Write (*, '(a)') ' ' // trim(s)
Stop
End If
!
! Number seems OK. Form it.
!
If (jtop==0) Then
rexp = 1.0D+00
Else
If (jbot==1) Then
rexp = 10.0D+00**(jsgn*jtop)
Else
rexp = 10.0D+00**(dble(jsgn*jtop)/dble(jbot))
End If
End If

r8 = dble(isgn)rexprtop/rbot

Return
End Subroutine s_to_r8
Subroutine timestamp()

!*********************************************************************72
!
!c TIMESTAMP prints out the current YMDHMS date as a timestamp.
!
! Parameters:
!
! None
!
Implicit None

Character *(8) ampm
Integer d
Character *(8) date
Integer h
Integer m
Integer mm
Character *(9) month(12)
Integer n
Integer s
Character *(10) time
Integer y

Save month

Data month/'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', 'September', &
'October ', 'November ', 'December '/

Call date_and_time(date, time)

Read (date, '(i4,i2,i2)') y, m, d
Read (time, '(i2,i2,i2,1x,i3)') h, n, s, mm

If (h<12) Then
ampm = 'AM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Noon'
Else
ampm = 'PM'
End If
Else
h = h - 12
If (h<12) Then
ampm = 'PM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Midnight'
Else
ampm = 'AM'
End If
End If
End If

Write (*, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d,&
month(m), y, h, ':', n, ':', s, '.', mm, ampm

Return
End Subroutine timestamp
!***************************************************************
!***************************************************************
!c $$$ Run Output:
!***************************************************************
!***************************************************************
set term png
set grid
set style data lines
unset key
set output "eot.png"
set xlabel "<---Normalized Date--->"
set ylabel "<---Minutes Early/Late--->"
set "The equation of time"
plot "analemma_data.txt" using 1:2 with lines
set output "declination.png"
set xlabel "<---Normalized Date--->"
set ylabel "<---Degrees North/South--->"
set "Declination"
plot "analemma_data.txt" using 1:3 with lines
set output "analemma.png"
set xlabel "<---Minutes Early/Late--->"
set ylabel "<---Degrees North/South--->"
set "The analemma"
plot "analemma_data.txt" using 2:3 with lines
quit

0.000000000 -4.42835237 -22.8216443
0.393316000E-04 -4.43494921 -22.8202411
0.786632001E-04 -4.44154475 -22.8188363
0.117994800E-03 -4.44813899 -22.8174299
0.157326400E-03 -4.45473192 -22.8160220
..................................
0.149857286E-01 -6.83471800 -22.1754863
0.150250633E-01 -6.84071876 -22.1735002
0.150643979E-01 -6.84671771 -22.1715126
0.151037327E-01 -6.85271484 -22.1695235
0.151430674E-01 -6.85871016 -22.1675329
...................................
0.799345326 15.4248445 -11.0310938
0.799385133 15.4270439 -11.0362340
0.799424939 15.4292410 -11.0413737
0.799464744 15.4314356 -11.0465128
0.799504550 15.4336279 -11.0516513
................................
0.999646016 -4.36892248 -22.8342033
0.999685347 -4.37553097 -22.8328141
0.999724679 -4.38213817 -22.8314233
0.999764010 -4.38874407 -22.8300310
0.999803342 -4.39534869 -22.8286371
0.999842674 -4.40195202 -22.8272417
0.999882005 -4.40855405 -22.8258447
0.999921337 -4.41515478 -22.8244461
0.999960668 -4.42175422 -22.8230460

ANALEMMA
FORTRAN95,2013&2018 version
Compute and plot the analemma, equation of time,
and declination.
This program is d on a C program by Brian Tung.

Created data file "analemma_data.txt".
Created command file "analemma_commands.txt".

ANALEMMA
Normal end of execution.
6 April 2022
...Program finished with exit code
Press ENTER to exit console.

地球仪8字曲线估算时间方程式公式差于一致的24小时和太阳的实际位置,
形成数据用于图形, 点数可上至五十万五千（505000）.运行验证数值演示如上，
并获得其图形命令，如果经处理应该可得示意表达.

广州