pro     makechiarr4_help
        print, 'makechiarr4,<file>,chi,e_rad,axratio,angle,corerad'
        print, 'marginalise4_1_2,chi,res'
end
 
pro     marginalise4_3_4,chi,res
        s = size(chi)
        sb = s[3]
        se = s[4]
        res = fltarr(sb,se)
        for i=0,sb-1 do begin
                for j = 0,se-1 do begin
                        res[i,j] = min(chi[*,*,i,j])
                endfor
        endfor
end

pro     marginalise4_1_2,chi,res
        s = size(chi)
        sb = s[1]
        se = s[2]
        res = fltarr(sb,se)
        for i=0,sb-1 do begin
                for j = 0,se-1 do begin
                        res[i,j] = min(chi[i,j,*,*])
                endfor
        endfor
end

; load the chi-squared data from filtered lensview log
; grep out lines with "main: Starting minimis" and "main: ImgDoConjGradient returned"
; from the log file. Then load them with this.
; Works for a single value of crit radius. You should chop up the filtered
;	file to unique vel dispersions
;	surf: the chi-squared value
; 	ordsh: the shear ordinate
;	ordsa: the shear angle ordinate
pro		makechiarr4, filename, surf, ordcritrad, ordsh, ordsa, ordcr
	restore,'/home/rwayth/idl/surfacetempl.dat'
	chidata = read_ascii(filename, TEMPLATE = template4)
; field 7	"lensoffset"
; field 8	"x:"
; field 9	x offset
; field 10	"offsety:"
; field 11	y offset
; field 12	"Critrad:" or some mass scale
; field 13	mass_scale
; field 14	"shear/Ellipticity"
; field 15	shear
; field 16	"shear_angle/Orient_Angle"
; field 17	shear angle (degrees)
; field 18	"core_rad"
; field 19	core radius
; field 20	"date"
; field 21	"time"
; field 22	"main:"
; field 23	"img_DoConjGradient"
; field 24	"returned"
; field 25	return value
; field 26	"Chi"
; field 27	"Squ:"
; field 28	chi-square value
; field 29	"Entropy:"
; field 30	entropy
; field 31,32 "Reduced chisqu:"
; field 33	reduced_chisqu
	xosize = size([uniq(chidata.FIELD09[sort(chidata.FIELD09)])])
	yosize = size([uniq(chidata.FIELD11[sort(chidata.FIELD11)])])
	mssize = size([uniq(chidata.FIELD13[sort(chidata.FIELD13)])])
	shsize = size([uniq(chidata.FIELD15[sort(chidata.FIELD15)])])
	sasize = size([uniq(chidata.FIELD17[sort(chidata.FIELD17)])])
	crsize = size([uniq(chidata.FIELD19[sort(chidata.FIELD19)])])
	shord = fltarr(shsize[1])
	saord = fltarr(sasize[1])
	critradord = fltarr(mssize[1])
	crord = fltarr(crsize[1])
print,'Mass scale size: ',strtrim(string(mssize[1]),2),', Shear size: ',strtrim(string(shsize[1]),2),'. Shear ang size: ',strtrim(string(sasize[1]),2),'. core rad size: ',strtrim(string(crsize[1]),2),'.'
	surface = fltarr(mssize[1], shsize[1],sasize[1],crsize[1])
	for k=0, mssize[1]-1 do begin
		critradord[k] = chidata.FIELD13[k*sasize[1]*shsize[1]*crsize[1]]
		for l=0, shsize[1] -1 do begin
			shord[l] = chidata.FIELD15[l*sasize[1]*crsize[1]]
			for i=0,sasize[1]-1 do begin
				saord[i] = chidata.FIELD17[i*crsize[1]]
				for j=0, crsize[1]-1 do begin
					crord[j] = chidata.FIELD19[j]
					surface[k,l,i,j] = chidata.FIELD28[k*sasize[1]*shsize[1]*crsize[1] + l*crsize[1]*sasize[1] + i*crsize[1] + j]	
				endfor
			endfor
		endfor
	endfor
	surf = surface
	ordsh = shord
	ordsa = saord
	ordcritrad = critradord
	ordcr = crord
end