;+ ; NAME: ; make_integral_counts ; ; PURPOSE: ; Compute the integral counts from the data "cube" ; ; CATEGORY: ; ; ; CALLING SEQUENCE: ; make_integral_counts ; ; INPUTS: ; wavelength ; save file IDL ; ; OUTPUTS: ; Flux (Jy) ; Integral number counts (/sr) ; Differential number counts (/sr) ; ; PROCEDURE: ; None ; ; EXAMPLE: ; make_integral_counts, 170.e-6, 1, Sjy_170, N_170, dnds_170, file= 'create_counts_0170_Omega_lambda.save', /ver ; ; MODIFICATION HISTORY: ; 12 feb-2000: written G. lagache ; 18 May 2002: make up for the delivery, G. Lagache ;- PRO make_integral_counts, wave, wnum, Sjy, N, dnds, file=file, plot =plot, STOP=STOP, restricted_z =restricted_z, z_in = z_in, verb=verb loadct_plot ; Restore arrays ;--------------- IF KEYWORD_SET(file) THEN BEGIN restore, file ENDIF ELSE BEGIN print, 'Please, enter your file name' ENDELSE IF KEYWORD_SET(restricted_z) AND KEYWORD_SET(z_in) THEN BEGIN good=where(z GT z_in - 0.15d0 *z and z LT z_in+0.15d0 *z) z=z(good) DNDLNLDZ = DNDLNLDZ(*, good) DSDZ = DSDZ(*, good) SLZ = SLZ(*, good) ENDIF ; COMPUTE DNDS ;-------------- delta_lnlum = (-alog(lum_array) + shift(alog(lum_array),-1))(0) dsdz_new = dsdz For il =0, n_elements(lum_array)-1 DO dsdz_new(il, *) = smooth(abs(dsdz(il,*)), 15) ; smooth discontinuities dndlnlds = dndlnldz / dsdz_new/ 1.d0 bad=where(dsdz EQ 0 OR finite(dndlnldz) EQ 0, count_bad) IF count_bad GT 0 THEN dndlnlds(bad)=0.d0 ; Creation of flux table ;----------------------- s_start = -7. ; start in log(Jy) s_end = 3. ; end in log(Jy) s_step = .05 ; Prendre un pas petit pour "IF (Sjy(i) GE S1(il) and sJy(i) LE S2(il))" ; mais pas en S doit etre plus grand que pas en L quand discontinuites s_loopend = FIX((s_end - s_start) / s_step) sjy = FLTARR(s_loopend +1) FOR is = 0, s_loopend DO BEGIN sjy(is) = s_start + is * s_step ENDFOR delta_s = sJy *0. for i=0, n_elements(sJy) -2 do delta_s(i) = - 10^sJy(i) + 10^sJy(i+1) ; FIND POINTS WHERE THE DS/DZ MAY SHOW DISCONTINUITIES ; (Case for large wavelegnths) ;----------------------------------------------------- z1 = fltarr(n_elements(lum_array)) z2 = z1 S2 = z1 S1 = z1 FOR ilum = 0, n_elements(lum_array) -1 DO BEGIN toto=reform(dsdz(ilum, *)) tutu=toto For i=0, n_elements(z)-1 DO tutu(i) = abs(toto(i)) -toto(i) g=where(tutu NE 0, count) IF count GT 1 THEN BEGIN ; Si 1 seul point negatif: c'est ignore z1(ilum) = z(g(0)) z2(ilum) = z(g(count-1)) S1(ilum) = alog10(Slz(ilum, g(0))) S2(ilum) = alog10(Slz(ilum, g(count-1))) ENDIF ENDFOR ; COMPUTE dN/dS's ;-------------------- dnds=sjy*0 & dn1 =0 & dn2 =0 & dn3 =0 FOR i =0, n_elements(sJy) -2 do BEGIN ; LOOP on fluxes bv =0*1.d0 FOR il = 0, N_elements(lum_array) -1 do BEGIN ; LOOP on L IF z1(il) EQ 0 THEN BEGIN ; PAS DE CHANGEMENT DE PENTE DANS LE SLz(il, *) temp_s = alog10(slz(il, *)) good = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count) IF count GT 0 THEN bv = bv + dndlnlds(il,good)* delta_lnlum ENDIF ELSE BEGIN IF z1(il) EQ z2(il) THEN BEGIN ; 1 SLOPE CHANGE PRINT, 'I never encoutered this case. to be checked' ;(check that S2 > S1 in that case and z2 = max(z)) STOP IF Sjy(i) GT S2(il) THEN BEGIN temp_s = alog10(slz(il, *)) good = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count) IF count GT 0 THEN bv = bv + dndlnlds(il,good)* delta_lnlum ENDIF IF SJy(i) LE S2(il) THEN BEGIN good1 = where(z LT z1(il)) temp_s = alog10(slz(il, good1)) & dn1= dndlnlds(il,good1) good1 = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count1) IF count1 GT 0 THEN dn1= dn1(good1) ELSE dn1 =0 good2 = where(z GE z1(il)) temp_s = alog10(slz(il, good2)) & dn2= dndlnlds(il,good2) good2 = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count2) IF count2 GT 0 THEN dn2= dn2(good2) ELSE dn2 =0 bv = bv + (dn1 + dn2)* delta_lnlum ENDIF ENDIF ELSE BEGIN IF z1(il) LT z2(il) AND z2(il) LE max(z) THEN BEGIN ; 2 SLOPE CHANGES IF (Sjy(i) GT S2(il) OR Sjy(i) LT S1(il)) THEN BEGIN ; in the region of S where slope is OK temp_s = alog10(slz(il, *)) ; get fluxes for this luminosity for all z good = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count) IF count GT 0 THEN bv = bv + dndlnlds(il,good)* delta_lnlum ENDIF IF (Sjy(i) GE S1(il) and sJy(i) LE S2(il)) THEN BEGIN ; in the region of S where slope change sign good1 = where(z LT z1(il)) ; find z below change of slope in z temp_s = alog10(slz(il, good1)) & dn1= dndlnlds(il,good1) ; get fluxes for this luminosity for z below good1 = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count1) IF count1 GT 0 THEN dn1= dn1(good1) ELSE dn1 =0 good2 = where(z GE z1(il) and z LT z2(il)) temp_s = alog10(slz(il, good2)) & dn2= dndlnlds(il,good2) ; find z in the region of change of slope in z good2 = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count2) IF count2 GT 0 THEN dn2= dn2(good2) ELSE dn2 =0 good3 = where(z GE z2(il)) ; find z after change of slope in z temp_s = alog10(slz(il, good3)) & dn3= dndlnlds(il,good3) good3 = where(abs(sJy(i)-temp_s) Eq min(abs(sJy(i)-temp_s)) AND Sjy(i) GE min(temp_s) AND Sjy(i) LE max(temp_s), count3) IF count3 GT 0 THEN dn3= dn3(good3) ELSE dn3 =0 bv = bv + (dn1 + dn2 + dn3)* delta_lnlum ENDIF ENDIF ELSE BEGIN print, 'Error' ENDELSE ENDELSE ENDELSE ENDFOR dnds(i) = bv ; = dn/ds ENDFOR ; Integral of the counts ;----------------------- N= fltarr(n_elements(sJy)) FOR i=0, n_elements(sJy) -2 DO BEGIN good=where(sJy GT sJy(i)) N(i)= total(dnds(good)*delta_s(good)) ENDFOR ; PLOT ;----- IF KEYWORD_SET(plot) THEN BEGIN title='' window, wnum, xs=400, ys=400 !p.charsize=1.2 plot, sJy, alog10(N), xtitle='Log (S) Jy', ytitle='N(>S)', title=title, xr=[s_start,s_end], /xsty oplot, sJy, alog10(N), color=1, thick=2 ENDIF ; CIB ;---- fond=total(dnds*delta_s*10.^sjy) * 3.e8/wave * 1.e-26 ; en W/m2/sr IF KEYWORD_SET(verb) THEN BEGIN print, 'CFIRB', fond, ' W/m2/sr' print, 'CFIRB', total(dnds*delta_s*10.^sjy)/1.e6, ' MJy/sr' ENDIF IF KEYWORD_SET(STOP) THEN STOP END