calculate equivalent width using python code

Question

I have this Fortran program for compute equivalent width of spectral lines i hope to find help for write python code to do same algorithm (input file contain tow column wavelength and flux)

     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

Ha! You asked for us to do it for you essentially... Provide your attempt and we'll shoot you some tips — FirebladeDan, Dec 03 '15 at 20:05

martineau · Answer 1 · 2015-12-20T03:57:45.293

Here's a fairly literal translation:

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()

However this would be a more idiomatic translation:

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()

score 0 · Answer 2 · edited Sep 21 '20 at 21:56

I would suggest that a more natural way to do this in Python is to focus on the properties of the spectrum itself, and use your parameters in astropy's specutils.

In particular equivalent_width details are here. For more general info on specutils, specutils.analysis and its packages follow these links:

specutils top level and specutils.analysis

To use this package you need to create a Spectrum1D object, the first component of which will be your wavelength axis and the second will be the flux. You can find details of how to create a Spectrum1D object by following the link in the analysis page (at the end of the third line of first paragraph).

It's a very powerful approach and has been developed by astronomers for astronomers.

calculate equivalent width using python code

2 Answers2