function arithtogram,x1,x2,xout,oper,plus=plus,minus=minus,divide=divide,\$ w1=w1,w2=w2,nbin=nbin,xmin=xmin,xmax=xmax,nonorm=nonorm, _extra=e ;+ ;function arithtogram ; return the result of an arithmetical operation on the ; frequency distributions of two lists ; ;syntax ; hx=arithtogram(x1,x2,xout,operator,/plus,/minus,/divide,\$ ; w1=w1,w2=w2,nbin=nbin,xmin=xmin,xmax=xmax,/nonorm) ; ;parameters ; x1 [INPUT; required] list of numbers whose frequency ; distribution must be multiplied with that of X2 ; x2 [INPUT; required] list of numbers whose frequency ; distribution must be multiplied with that of X1 ; * it is assumed that X1 and X2 both have the same units ; because otherwise this kind of multiplication makes ; no sense ; xout [OUTPUT; required] the list of numbers at which the ; output is tabulated ; * by default, this is the sorted unique set of [X1,X2] ; * if NBIN is set, generates a logarithmic or linear ; grid of that size ; oper [INPUT] must be one of '+', '-', '/', '*': ; the arithmetical operation to carry out on the ; frequency distributions of X1 and X2, i.e., ; result = X1 OPER X2 ; * if not given, assumed to be '*', unless one of ; keywords PLUS, MINUS, or DIVIDE are set ; ;keywords ; plus [INPUT] if set, OPER is assumed to be '+' ; minus [INPUT] if set, OPER is assumed to be '-' ; divide [INPUT] if set, OPER is assumed to be '/' ; * DIVIDE takes precedence over MINUS takes precedence ; over PLUS ; w1 [INPUT] optional weight to be applied to histogram(X1) ; w2 [INPUT] optional weight to be applied to histogram(X2) ; * if not given, W1 and W2 are assumed to be 1 ; nbin [INPUT] number of bins in the output ; * if not set, will be the sorted unique set of [X1,X2] ; * if -ve, will produce logarithmic gridding ; xmin [INPUT] minimum value in the output grid ; * by default, uses min(X1)>min(X2) ; xmax [INPUT] maximum value in the output grid ; * by default, uses max(X1) min(X2) xmax=max(X1) < max(X2) mbin=0L & if keyword_set(nbin) then mbin=long(nbin[0]) wt1=1. & wt2=1. if keyword_set(w1) then wt1=w1[0] if keyword_set(w2) then wt2=w2[0] ; output grid if keyword_set(mbin) then begin if mbin lt 0 then begin lxmin=alog10(xmin) & lxmax=alog10(xmax) dx=(lxmax-lxmin)/(abs(mbin)-1L) xout=10.^(findgen(abs(mbin)+1L)*dx+lxmin) endif else begin dx=(xmax-xmin)/(mbin-1L) xout=findgen(mbin+1L)*dx+xmin endelse endif else begin xout=[X1,X2] & xout=xout[uniq(xout,sort(xout))] ok=where(xout ge xmin and xout le xmax,mok) if mok eq 0 then message,'BUG!' xout=xout[ok] endelse ; make cdfs o1=sort(X1) & o2=sort(X2) c1=findgen(n1)/(n1-1L) & c2=findgen(n2)/(n2-1L) cc1=interpol(c1,x1[o1],xout) & cc2=interpol(c2,x2[o2],xout) dc1=n1*(cc1[1:*]-cc1) & dc2=n2*(cc2[1:*]-cc2) dxout=xout[1:*]-xout & xout=0.5*(xout[1:*]+xout) if not keyword_set(nonorm) then begin dc1=dc1/dxout & dc2=dc2/dxout endif case op of '+': hx=(wt1*dc1)+(wt2*dc2) '-': hx=(wt1*dc1)-(wt2*dc2) '/': hx=(wt1*dc1)/(wt2*dc2) else: hx=(wt1*dc1)*(wt2*dc2) endcase return,hx end