-5

我有这个 Fortran 程序来计算谱线的等效宽度我希望找到帮助编写 python 代码来执行相同的算法（输入文件包含两列波长和通量） 

     PARAMETER (N=195)              ! N is the number of data
     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
     DIMENSION X(N),Y(N)
     OPEN(1,FILE='halpha.dat')
     DO 10 I=1,N
     READ(1,*)X(I),Y(I)
     WRITE(*,*)X(I),Y(I)
10   CONTINUE
     CALL WIDTH(X,Y,N,SUM)     
     WRITE(*,*)SUM
     END
 c-----------------------------------------
    SUBROUTINE WIDTH(X,Y,N,SUM)
    PARAMETER (NBOD=20000)
    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    DIMENSION X(NBOD),Y(NBOD) 
    SUM=0.D0
    DO I=2,N
    SUM=SUM+(X(I-1)-X(I))*((1.-Y(I-1))+(1.-Y(I)))
C   WRITE(*,*)SUM
    END DO 
    SUM=0.5*dabs(SUM)
    RETURN
    END
4

2 回答 2

2

这是一个相当直译的翻译：

def main():
    N = 195  # number of data pairs
    x, y = [0 for i in xrange(N)], [0 for i in xrange(N)]

    with open('halpha.dat') as f:
        for i in xrange(N):
            x[i], y[i] = map(float, f.readline().split())
            print x[i], y[i]

    sum = width(x, y, N)
    print sum

def width(x, y, N):
    sum = 0.0
    for i in xrange(1, N):
        sum = sum + (x[i-1] - x[i]) * ((1. - y[i-1]) + (1. - y[i]))
    sum = 0.5*abs(sum)
    return sum

main()

然而，这将是一个更惯用的翻译：

from math import fsum  # more accurate floating point sum of a series of terms

def main():
    with open('halpha.dat') as f:  # Read file into a list of tuples.
        pairs = [tuple(float(word) for word in line.split()) for line in f]

    for pair in pairs:
        print('{}, {}'.format(*pair))

    print('{}'.format(width(pairs)))

def width(pairs):
    def term(prev, curr):
        return (prev[0] - curr[0]) * ((1. - prev[1]) + (1. - curr[1]))

    return 0.5 * abs(fsum(term(*pairs[i-1:i+1]) for i in range(1, len(pairs))))

main()
于 2015-12-03T21:00:03.090 回答
0

我建议在 Python 中执行此操作的一种更自然的方法是关注光谱本身的属性，并在 astropy 的 specutils 中使用您的参数。

特别是等效宽度细节在这里。有关 specutils、specutils.analysis 及其包的更多一般信息，请点击以下链接：

specutils 顶级 和 specutils.analysis

要使用这个包，您需要创建一个 Spectrum1D 对象，其第一个组件将是您的波长轴，第二个将是通量。您可以通过分析页面中的链接（在第一段第三行的末尾）找到有关如何创建 Spectrum1D 对象的详细信息。

这是一种非常强大的方法，由天文学家为天文学家开发。

于 2020-09-21T17:20:34.953 回答