;Spectralmap is a program that takes a number of total intensity radio images,
;along with frequency and rms information about the maps, and creates a spectral
;index map and a chi^2 map.  The input text file FREQDATA should have the frequency 
;of each map, map name, and off source rms.  For example:
;  freq(Hz)         mapname       rms(Jy/beam)
;  350E+06          map1.fits     3.25E-04   
;  ...              ...           ...
;  ...              ...           ...
;Start IDL in the directory with the maps and this text file and type 
;>spectralmap
;and the maps 'SPECTRALINDEX.FITS' and 'CHISQRD.FITS' will be created inn that
;directory (it may take a few seconds).  I have it here for 8 maps, but I'll soon
;fix it for any number of maps (it is easy to do if you need to). 


pro spectralmap

readcol, 'FREQDATA', format= 'd,a,d', freq, iname, rms_vector
n=size(freq,/n_elements)  ;number of maps
lgfreq=alog10(freq)

spec_image=dblarr(1024,1024)
chi_image=dblarr(1024,1024)
spec_vector=dblarr(8)

i1=readfits(iname[0])
i2=readfits(iname[1])
i3=readfits(iname[2])
i4=readfits(iname[3])
i5=readfits(iname[4])
i6=readfits(iname[5])
i7=readfits(iname[6])
i8=readfits(iname[7])


for i=0, 1023 do begin
for j=0, 1023 do begin

temp1=i1[i,j]
temp2=i2[i,j]
temp3=i3[i,j]
temp4=i4[i,j]
temp5=i5[i,j]
temp6=i6[i,j]
temp7=i7[i,j]
temp8=i8[i,j]

if (temp1 ge 20*rms_vector[0]) and (temp2 ge 20*rms_vector[1]) and (temp3 ge 20*rms_vector[2]) and (temp4 ge 20*rms_vector[3]) and (temp5 ge 20*rms_vector[4]) and (temp6 ge 20*rms_vector[5]) and (temp7 ge 20*rms_vector[6]) and (temp8 ge 20*rms_vector[7]) then begin 

spec_vector[0]=alog10(temp1)
spec_vector[1]=alog10(temp2)
spec_vector[2]=alog10(temp3)
spec_vector[3]=alog10(temp4)
spec_vector[4]=alog10(temp5)
spec_vector[5]=alog10(temp6)
spec_vector[6]=alog10(temp7)
spec_vector[7]=alog10(temp8)

si=linfit(lgfreq,spec_vector,chisq=chi)

spec_image[i,j]=-si[1]
chi_image[i,j]=chi

endif else begin

spec_image[i,j]=0
chi_image[i,j]=0

endelse

endfor
endfor

writefits, 'SPECTRALINDEX.FITS', spec_image
writefits, 'CHISQRD.FITS', chi_image

end