;This program will take sets of Q and U maps, one set for each channel, and perform RM-synth on them.
;You need two input text files. 'Pfiles' is a three column text file with the Q fits file names, the
;U fits file names, and the central wavelength of the channels (I will soon change this part to deal
;in frequency space instead). 'Faradaydepth' is a two column text file that has the output Faraday
;dispersion fuction file names and the faraday depth (the Faraday dispersion fuction is the value of a
;given pixel as a fuction of Faraday depth).

pro rm_synthnew

readcol, 'filenames', format= 'a,a,d', qfiles, ufiles, lambda   ;read FITS file names and band wavelenth centers into arrays
readcol, 'FaradayDepth', format = 'a,d', fitsnames, fardepth   ;read range of faraday depths into an array
n=size(qfiles,/n_elements)  ;number of channels
b=size(fardepth,/n_elements)    ;number of faraday depths to try
lambdanot=mean(lambda)    ;the weighted average of the observed lambda^2  

qmaps=fltarr(1024,1024,n)
umaps=fltarr(1024,1024,n)

for i=0, n-1 do begin

qtemp = readfits(qfiles[i])
utemp = readfits(ufiles[i])
qmaps[*,*,i]=qtemp
umaps[*,*,i]=utemp

endfor

for j=0, b-1 do begin
f=0
for i=0, n-1 do begin

k=i

qtm=qmaps[*,*,k]
utm=umaps[*,*,k]

f = f + complex(qtm,utm)*complex(cos(-2*fardepth[j]*(lambda[k]-lambdanot)),sin(-2*fardepth[j]*(lambda[k]-lambdanot)))

endfor

final=sqrt(float(f*conj(f)))/n

writefits, fitsnames[j], final

endfor

end