I have trouble passing 2D arrays to fortran. I want to combine a bunch of not overlapping spectra. First I select the points on the x-axis, then I interpolate all data to this new, common grid. I store the spectra in a 2D list in python.
This works in Python 2.7, but very slow:
for i in range(len(wlp)):
print wlp[i],
for a in range(len(datax)):
inrange = 0
if datax[a][0] >= wlp[i] or datax[a][-1] <= wlp[i]:
for b in range(len(datax[a])-1):
if float(datax[a][b]) <= wlp[i] and float(datax[a][b+1]) >= wlp[i]:
sp = float(datax[a][b]); ep = float(datax[a][b+1])
delx = ep-sp; dely = float(data[a][b+1])-float(data[a][b])
ji = (dely/delx)*(wlp[i]-sp)+float(data[a][b])
inrange = 1
if inrange == 0: ji = '?0'
else: ji = ji * weights[a]
print ji,
print
The common x-grid is printed in column one and all the interpolated spectra are printed in subsequent columns. If there are some shorter ones out of range, it prints "?0". This helps to set up proper weights for each datapoints later.
I ended up having this fortran subroutine to speed it up with f2py:
c wlp = x axis points (wavelength)
c lenwlp = length of list wlp, len(wlp)
c datay = 2D python list with flux
c datax = 2D python list with wavelength
c lendatax = number of spectra, len(datax)
c datax_pl = list of the lengths of all spectra
c weights = list of optional weights
c maxi = length of the longest spectrum
C============================================================================80
SUBROUTINE DOIT(wlp,lenwlp,datay,datax,lendatax,datax_pl,
. weights,maxi)
C============================================================================80
INTEGER I,a,b,lenwlp,inrange,datax_pl(*),maxi,lendatax
DOUBLE PRECISION WLP(*),SP,EP,DELY,DELX,ji
DOUBLE PRECISION WEIGHTS(*)
DOUBLE PRECISION DATAY(lendatax,maxi)
DOUBLE PRECISION DATAX(lendatax,maxi)
2 FORMAT (E20.12, 1X, $)
3 FORMAT (A, $)
4 FORMAT (1X)
I = 1
DO WHILE (I.LE.lenwlp)
WRITE(*,2) WLP(I)
DO a=1,lendatax
inrange = 0
ji = 0.0
IF (datax(a,1).ge.WLP(I) .or.
. datax(a,datax_pl(a)).le.WLP(I)) THEN
DO b=1,datax_pl(a)-1
IF (DATAX(a,b).LE.WLP(I) .and.
. DATAX(a,b+1).GE.WLP(I)) THEN
SP = DATAX(a,b); EP = DATAX(a,b+1)
DELX = EP - SP; DELY = datay(a,b+1)-datay(a,b)
if (delx.eq.0.0) then
ji = datay(a,b)
else
ji = (DELY/DELX)*(WLP(I)-SP)+datay(a,b)
end if
inrange = 1
END IF
END DO
END IF
IF (inrange.eq.0) THEN
WRITE(*,3) ' ?0'
ELSE
WRITE(*,2) ji*WEIGHTS(a)
END IF
END DO
I = I + 1
write(*,4)
END DO
END
which compiles with gfortran 4.8 fine. Then I import it in the Python code, set up the lists and run the subroutine:
import subroutines
wlp = [...]
data = [[...],[...],[...]]
datax = [[...],[...],[...]]
datax_pl = [...]
weights = [...]
maxi = max(datax_pl)
subroutines.doit(wlp,len(wlp),data,datax,len(datax),datax_pl,weights,maxi)
and it returns:
ValueError: setting an array element with a sequence.
I pass the lists and the length of the longest spectrum (maxi), this should define the maximum dimension in fortran (?). I don't need return values, everything is printed on stdout.
The problem must be right at the beginning at the array declarations. I don't have experience with this... any advice is appreciated.
