;***************************************************************************** ;* program from Jack Scudder - for detector intercallibration work. ;***************************************************************************** PRO find_con_errors, c68, c95, c99, f1new, sigx, sigy, x, f2new, y, $ nparms, iswitch, status, rchi2 ;calculate the curvature matrix exactly ; Renormalize the matrix if rchi2 < 1 - gives better error estimates IF (rchi2 LT 1.) THEN renorm = rchi2 ELSE renorm = 1. denom=(f1new*sigx)^2+sigy^2 num=f1new*x+f2new-y m11=0.5d0/renorm* total( 2*x^2/(denom) - $ 2*(num)*x*2*f1new* sigx^2/(denom)^2 - $ (2*(sigx*num)^2/(denom)^2)- $ (2*f1new*sigx^2*2*(num)*x/(denom)^2) + $ 2*2*f1new*(sigx*num)^2*2*f1new*sigx^2/(denom)^3 ) m12=0.5d0/renorm* total ( 2*x/(denom) - $ 2*2*f1new*(num)*1*(sigx/denom)^2 ) m22=0.5d0/renorm* total( 2/(denom)) ; calculate the approximate error boxes about the best fit IF (nparms EQ 1) THEN BEGIN thetap = [ 1.0d0, 3.84d0, 6.63d0 ] ENDIF ELSE BEGIN thetap = [ 2.3d0, 5.99d0, 9.21d0 ] ENDELSE IF ( iswitch eq 1 ) then begin IF (status ne 4) then begin errorf195=(m22*thetap/(m11*m22-m12^2))^0.5 errorf295=(m11*thetap/(m11*m22-m12^2))^0.5 c68 = [errorf195(0),errorf295(0)] ; errors [slope,intercept] c95 = [errorf195(1),errorf295(1)] c99 = [errorf195(2),errorf295(2)] endif else begin ; ireturn code =4 exit error295= (thetap/m22)^0.5 c68 = [0,error295(0)] c95 = [0,error295(1)] c99 = [0,error295(2)] endelse endif else begin IF (status ne 4) then begin errorf195=(thetap/m11)^0.5 c68 = [errorf195(0),0] c95 = [errorf195(1),0] c99 = [errorf195(2),0] endif else begin ; ireturn code =4 exit print, 'I''m here and I don''t know why' help, iswitch, status stop endelse ENDELSE return END ;------------------------------------------------------------- PRO find_ks_con_errors, ks_c68, ks_c95, ks_c99, nparms, f1, $ save1_kappa, kappa, kappa_data, nsize, kappa_picture, $ kappa_temp=kappa_temp, titstr=titstr ; save1_kappa is the minimum index of kappa. IF (keyword_set(kappa_temp)) THEN chi2 = kappa_temp ELSE chi2 = kappa ks_c68 = dblarr(2) ks_c95 = dblarr(2) ks_c99 = dblarr(2) ; Define ks_c68, ks_c95 and ks_c99 from the kappa search ; These error estimates may be asymmetric about the slope IF (nparms EQ 1) THEN con68 = 1.00 ELSE con68 = 2.30 IF (nparms EQ 1) THEN con95 = 3.84 ELSE con95 = 5.99 IF (nparms EQ 1) THEN con99 = 6.63 ELSE con99 = 9.21 ; Determine how 'strong' the slope is IF (f1(save1_kappa) GT 0) THEN BEGIN ; min on the right IF (chi2(0) LT chi2(save1_kappa)+con68) THEN $ kappa_data.both68 = 1 IF (chi2(0) LT chi2(save1_kappa)+con95) THEN $ kappa_data.both95 = 1 ENDIF ELSE BEGIN ; min on the left IF (chi2(nsize-1) LT chi2(save1_kappa)+con68) THEN $ kappa_data.both68 = 1 IF (chi2(nsize-1) LT chi2(save1_kappa)+con95) THEN $ kappa_data.both95 = 1 ENDELSE i68 = 0 i95 = 0 i99 = 0 ; march to the right to find where we cross the confidence level FOR ii=save1_kappa+1,n_elements(f1)-1 DO BEGIN IF (i68 EQ 1 AND i95 EQ 1 AND i99 EQ 1) THEN GOTO, out_now IF ( (chi2(ii)-chi2(save1_kappa)) GE con68 AND i68 EQ 0) $ THEN BEGIN alpha= $ ( ( chi2(save1_kappa)+ con68 ) - chi2(ii-1) ) / $ ( chi2(ii) - chi2(ii-1) ) intersect= alpha * (f1(ii)-f1(ii-1)) + f1(ii-1) ks_c68(1) = abs(intersect-f1(save1_kappa)) i68 = 1 ENDIF IF ( (chi2(ii)-chi2(save1_kappa)) GE con95 AND i95 EQ 0) $ THEN BEGIN alpha= $ ( ( chi2(save1_kappa)+ con95 ) - chi2(ii-1) ) / $ ( chi2(ii) - chi2(ii-1) ) intersect= alpha * (f1(ii)-f1(ii-1)) + f1(ii-1) ks_c95(1) = abs(intersect-f1(save1_kappa)) i95 = 1 ENDIF IF ( (chi2(ii)-chi2(save1_kappa)) GE con99 AND i99 EQ 0) $ THEN BEGIN alpha= $ ( ( chi2(save1_kappa)+ con99 ) - chi2(ii-1) ) / $ ( chi2(ii) - chi2(ii-1) ) intersect= alpha * (f1(ii)-f1(ii-1)) + f1(ii-1) ks_c99(1) = abs(intersect-f1(save1_kappa)) i99 = 1 ENDIF ENDFOR out_now: ; march to the left to find where we cross the confidence level i68 = 0 i95 = 0 i99 = 0 FOR ii=save1_kappa-1,0,-1 DO BEGIN IF (i68 EQ 1 AND i95 EQ 1 AND i99 EQ 1) THEN GOTO, out_now2 IF ( (chi2(ii)-chi2(save1_kappa)) GE con68 AND i68 EQ 0) $ THEN BEGIN alpha= $ ( ( chi2(save1_kappa)+ con68 ) - chi2(ii+1) ) / $ ( chi2(ii) - chi2(ii+1) ) intersect= alpha * (f1(ii)-f1(ii+1)) + f1(ii+1) ks_c68(0) = abs(f1(save1_kappa)-intersect) i68 = 1 ENDIF IF ( (chi2(ii)-chi2(save1_kappa)) GE con95 AND i95 EQ 0) $ THEN BEGIN alpha= $ ( ( chi2(save1_kappa)+ con95 ) - chi2(ii+1) ) / $ ( chi2(ii) - chi2(ii+1) ) intersect= alpha * (f1(ii)-f1(ii+1)) + f1(ii+1) ks_c95(0) = abs(f1(save1_kappa)-intersect) i95 = 1 ENDIF IF ( (chi2(ii)-chi2(save1_kappa)) GE con99 AND i99 EQ 0) $ THEN BEGIN alpha= $ ( ( chi2(save1_kappa)+ con99 ) - chi2(ii+1) ) / $ ( chi2(ii) - chi2(ii+1) ) intersect= alpha * (f1(ii)-f1(ii+1)) + f1(ii+1) ks_c99(0) = abs(f1(save1_kappa)-intersect) i99 = 1 ENDIF ENDFOR out_now2: IF (kappa_picture) THEN BEGIN IF (keyword_set(kappa_temp)) THEN ytitle = 'kappa_temp' ELSE ytitle = 'kappa' xrange = [f1(0),f1(n_elements(f1)-1)] plot,f1,chi2,xtitle='f1',ytitle=ytitle,xstyle=1,ystyle=1,xrange=xrange,title=titstr,/ylog oplot, [f1(save1_kappa),f1(save1_kappa)],[min(chi2),max(chi2)],line=0 oplot, [f1(save1_kappa)-ks_c68(0),f1(save1_kappa)-ks_c68(0)],$ [min(chi2),max(chi2)],line=2 oplot, [f1(save1_kappa)+ks_c68(1),f1(save1_kappa)+ks_c68(1)],$ [min(chi2),max(chi2)],line=2 oplot, [f1(save1_kappa)-ks_c95(0),f1(save1_kappa)-ks_c95(0)],$ [min(chi2),max(chi2)],line=4 oplot, [f1(save1_kappa)+ks_c95(1),f1(save1_kappa)+ks_c95(1)],$ [min(chi2),max(chi2)],line=4 time_stamp ENDIF return END ; -------------------------------------------------------------- FUNCTION f2func,f1dum,x,y,sigx,sigy,iswitch f1dum = double(f1dum) x = double(x) y = double(y) sigx = double(sigx) sigy = double(sigy) answer = double(0.d0) answer = (-total((f1dum*x-y)/((f1dum*sigx)^2+sigy^2),/double)/$ total(1.0d0/((f1dum*sigx)^2+sigy^2),/double))*iswitch ;;ok jds return,answer END;========================================================= FUNCTION dkappadf1,f1dum,f2dum,x,y,sigx,sigy,iswitch f1dum = double(f1dum) f2dum = double(f2dum) x = double(x) y = double(y) sigx = double(sigx) sigy = double(sigy) denom = double(0.d0) testdum = double(0.d0) df2df1 = double(0.d0) dchidf2 = double(0.d0) answer = double(0.d0) ; print,format='(3(e25.18,2x))', x(0),x(1),x(2) ; print,format='(3(e25.18,2x))', y(0),y(1),y(2) ; print,format='(3(e25.18,2x))', sigx(0),sigx(1),sigx(2) ; print,format='(3(e25.18,2x))', sigy(0),sigy(1),sigy(2) ; print,format='(2(e25.18,2x))', f1dum,f2dum denom=(f1dum*sigx)^2+sigy^2 testdum= total (2.d0*(f1dum*x+f2dum-y)*x/(denom) - $ (f1dum*x+f2dum-y)^2*2.d0*f1dum*sigx^2/(denom)^2,/double );;; df2df1= - ( ($ total(1/denom,/double)*$ total(x/denom-2.d0*f1dum*sigx^2*(f1dum*x-y)/denom^2,/double)$ +total((f1dum*x-y)/denom,/double)* $ total(2.d0*f1dum*sigx^2/denom^2,/double) $ )/$ (total(1./denom,/double))^2 ) dchidf2 =total(2.d0*(f1dum*x+f2dum-y)/denom,/double) answer = testdum+ dchidf2*df2df1*iswitch ; print,format='(e25.18,2x)', answer return,answer END;===================================================================== ; -----DOCUMENTATION----- ;+ ;======================================================================== ; NAME: gen_fit_release.pro ; ; Purpose: To develop general fitting routine by ; hypothesis testing for straight line model ; with smart start and Marquardt back end ; including error estimates. ; ; Written, 5/1/97, Jack Scudder ; Modified, 5/2/97, Pamela Puhl-Quinn and Jeremy Faden ; 5/19/97, added guess_slope keyword ; 6/27/97, added intercept-only code JDS/RDH (gen_fit_mod) ; 8/11/97 placed in common directory and renamed gen_fit_v04 ; 1/17/98 renamed gen_fit_v05, fixed c95 dependence on ; nparms ; 5/14/98 - * Renamed gen_fit_release ; * Added force_kappa_search keyword: ; Kappa search method will be used to ; bracket the zero; error estimates in SLOPE ; ONLY will be made using the kappa ; surface itself; these estimates are ; returned in ks_68, ks_95 and ks_99. ; * Added n_iterations keyword ; * Cleaned up user interface (doc_library friendly) ; Just type gen_fit_release to see the documentation ; 10/15/98 - Refined the logic regarding determining errors ; when reduced chisq is < 1. Renormalize the chisq ; surface so that the minimum value is on the ; order of the # of degrees of freedom, THEN ; find errors on the fit parameters from this temporary ; chisq surface. This logic applies to the ; determination of c68, c95, c99, ks_68, ks_95, ; ks_99. ; 04/13/99 - Fixed the rchisq normalization - Wasn't correct ; for nparms not equal to 2. ; ; Curator, Pamela Puhl-Quinn, ppq@space-theory.physics.uiowa.edu ; ; ------------------- KEYWORDS ---------------------- ; ; INFO: Calls doc_library to display this header documentation. ; 1: print to terminal, 2: print to default printer ; ; x=x,sigma_x=sigma_x - abscissa data and errors (IN) ; y=y,sigma_y=sigma_y - ordinate data and errors (IN) ; titstr=titstr - title string for kappa picture (IN) ; n_iterations=n_iterations - Maximum # of iterations to use when ; Bolzano is zeroing in (default=100) (IN) ; nparms=nparms - number of parameters to fit (IN) ; 0 = intercept only (slope=1) ; 1 = slope only (intercept=0) ; 2 = slope and intercept (default) ; guess_slope=guess_slope - your best guess on the slope value (IN,optional) ; errorpercent=errorpercent - converge to this rel. tolerence (IN) ; .06 = agreement w/i 3 sig. figs ; .002 = agreement w/i 4 sig. figs (default) ; .0002 = agreement w/i 5 sig. figs ; .00004 = agreement w/i 6 sig. figs ; .000001= agreement w/i 7 sig. figs ; slope=slope - slope (OUT) ; intercept=intercept - intercept (OUT) ; chisq=chisq - chi-squared of the fit (OUT) ; rchisq=rchisq - reduced chi-squared of the fit (OUT) ; s2=s2 - (normalized distance of each point from fit) ^ 2 ; ; c68=c68, c95=c95, c99=c99 - 68%, 95% and 99% confidence ; level symmetric errors in slope AND intercept ; (OUT) ; e.g., c68=[slope error, intercept error] ; ; ks_c68=ks_c68,ks_c95=ks_c95,ks_c99=ks_c99 ; - 68%, 95% and 99% confidence ; level asymmetric error in SLOPE ONLY ; (OUT only when force_kappa_search is set) ; e.g., ks_c68=[slope error to left, slope error to right] ; NOTE: 68, 95 AND 99 CAN ALL BE IDENTICAL IF ; THE NUMBER OF DATA POINTS IS LARGE AND THE ; CHI-SQUARED IS NARROW (RUN INTO CHI-SQUARED ; GRID RESOLUTION LIMITATIONS). USE THE ; kappa_picture KEYWORD IF YOU WANT TO SEE THIS.. ; ; ks_range_deg=[x,y] specify range for kappa search. ; ; /force_kappa_search - Forces kappa search method to bracket zero; ; ks_c68, ks_c95, ks_c99 are filled with ; asymmetric slope errors determined from ; the kappa surface itself ; ; /verbose - Diagnostics written out ; ; kappa_picture=nwindow - window number for kappa vs f1 plot (nonzero) (IN) ; /print - In case you want to print the above plot ; itime=itime - Pam's diag parameter (IN) ; kappa_data=kappa_data - Data to make the kappa vs f1 picture (OUT) ; status=status - 2=fail, kappa search failed (OUT) ; 1=pass, zero bracketed,passed tolerence test on f1 ; 0=fail, dkappadf1 < 1.e-19 ; 3=fail, zero brack., passed tol., MAXIMUM! ; 4=pass, For nparms=0, slope=1 used, intercept found ; 5=fail, Couldn't converge within n_iterations ;======================================================================== ;- ; -----DOCUMENTATION----- PRO gen_fit_release, $ x=x,sigma_x=sigx,y=y,sigma_y=sigy, $ n_iterations=n_iterations, $ titstr=titstr, $ nparms=nparms, $ guess_slope=guess_slope,$ errorpercent=errorpercent, $ slope=slope, intercept=intercept, $ chisq=chisq, rchisq=rchisq, $ s2=s2, $ c68=c68, c95=c95, c99=c99,$ ks_c68=ks_c68, ks_c95=ks_c95, ks_c99=ks_c99,$ force_kappa_search=force_kappa_search,$ ks_range_degrees= ks_range_deg, $ renormalize=renormalize, $ verbose=verbose,$ kappa_picture=kappa_picture,$ print=print,$ itime=itime,$ kappa_data=kappa_data,$ status=status, $ info = info ;------------------------------------------------------------------- if n_elements( ks_range_deg ) eq 0 then begin ks_range_deg= [ -85, 85 ] endif IF keyword_set(INFO) THEN BEGIN IF info EQ 2 THEN print = 1 ELSE print = 0 doc_library, 'gen_fit_release', print = print return ENDIF IF (NOT keyword_set(x)) THEN BEGIN doc_library, 'gen_fit_release' return ENDIF IF (keyword_set(force_kappa_search)) THEN force_kappa_search = 1 ELSE force_kappa_search = 0 IF (keyword_set(kappa_picture)) THEN BEGIN n_picture = kappa_picture kappa_picture = 1 ENDIF ELSE BEGIN kappa_picture = 0 ENDELSE IF (keyword_set(verbose)) THEN verbose = 1 ELSE verbose = 0 IF (keyword_set(n_iterations)) THEN n_iterations = n_iterations ELSE n_iterations = 100 IF (keyword_set(titstr)) THEN titstr = titstr ELSE titstr = 'KAPPA PICTURE' IF (keyword_set(force_kappa_search)) THEN BEGIN ks_required = 1 IF (verbose OR kappa_picture) THEN print, 'Bracketing zero will be determined from a kappa search because you set force_kappa_search' ENDIF ELSE BEGIN ks_required = 0 ENDELSE IF (keyword_set(print)) THEN print = 1 ELSE print = 0 IF (print) THEN set_ps, orient='L' if n_elements(nparms) eq 0 then nparms = 2 if nparms lt 0 or nparms gt 2 then begin ;jds 6/27/97 print, 'one or two parameter fit only: nparms=1 or nparms=2' return endif if not keyword_set(errorpercent) then errorpercent=0.002d0 IF (nparms eq 2 OR nparms eq 0) then iswitch=1 ; jds 6/27/97 IF (nparms eq 1) then iswitch=0 IF (ks_required AND nparms EQ 0) THEN message, '/force_kappa_search does not work with nparms = 0' x=double(x) y=double(y) sigx=double(sigx) sigy=double(sigy) ireturncode = 99 zero=1d-19 npts= n_elements(x) ; nparms=0 means fit to a slope=1, find the best intercept IF (nparms EQ 0) THEN BEGIN ; jds 6/27/97 f1new = 1.0d0 f2new = f2func(f1new,x,y,sigx,sigy,iswitch) ireturncode = 4 goto, terminal ENDIF IF (keyword_set(guess_slope)) then begin slope1 = guess_slope intercept1 = f2func(slope1,x,y,sigx,sigy,iswitch) slope_init=slope1 intercept_init=intercept1 ENDIF ELSE BEGIN ; Assuming no error in x a=total((x/sigy)^2) b=total((1./sigy)^2) c=total(x/sigy^2) det=a*b-c^2 y1=total(x*y/sigy^2) y2=total(y/sigy^2) slope1= (b*y1-c*y2)/det intercept1= (-c*y1+a*y2)/det slope_init = slope1 intercept_init = intercept1 ENDELSE ; Now, calculate the dkappadf1 here at slope1 IF (verbose OR kappa_picture) THEN print, 'Initializing...' init_counter = 0 init_loop: f1dum = slope1 init_counter = init_counter + 1 ; Count attempts to bracket zero f2dum = f2func(f1dum,x,y,sigx,sigy,iswitch) testinit1 = dkappadf1(f1dum,f2dum,x,y,sigx,sigy,iswitch) intercept1 = f2dum ; Find slope2 from this derivative information kappaslope1 = total((f1dum*x + f2dum - y)^2/((f1dum*sigx)^2+sigy^2)) slope2 = slope1 - kappaslope1/testinit1 ; Now what is dkappadf1 at slope2? f1dum = slope2 f2dum = f2func(f1dum,x,y,sigx,sigy,iswitch) testinit2 = dkappadf1(f1dum,f2dum,x,y,sigx,sigy,iswitch) intercept2 = f2dum IF (verbose OR kappa_picture) THEN BEGIN print, 'Init_loop: slope1,slope2',slope1,slope2 print, 'Init_loop: intercept1,intercept2',intercept1,intercept2 print, 'Init_loop: testinit1,testinit2',testinit1,testinit2 ENDIF ; Bypass - If ks_required, we MUST do a kappa ; search. Therefore we're going to bypass the following logic about ; bracketing the zero and get right to the kappa search IF (ks_required) THEN GOTO, search_kappaspace ; Evaluate whether or not we've bracketed the zero: if (testinit1*testinit2 lt 0) then begin arrays = dblarr(2) arrayt = dblarr(2) arrayi = dblarr(2) arrays = [slope1,slope2] arrayt = [testinit1,testinit2] arrayi = [intercept1,intercept2] f1left = min(arrays,imin) f1right = max(arrays,imax) testleft = arrayt(imin) testright = arrayt(imax) f2left = arrayi(imin) f2right = arrayi(imax) IF (kappa_picture) THEN goto, plot_kappa goto, zero_in endif else begin ; dkappakf1 is virtually zero if(abs(testinit2) lt zero) then begin f1new = slope2 f2new = f2dum ireturncode = 0 GOTO, terminal endif ; Can't bracket the zero, too many attempts if (init_counter ge 1) then begin GOTO, search_kappaspace endif ; Try again to bracket the zero slope1=slope2 goto, init_loop ENDELSE ;===================KAPPA PICTURE (UNNECESSARY SEARCH)=========== plot_kappa: print, 'Performing kappa search for plotting only (unnecessary)' IF (NOT print) THEN window,n_picture,xsize=550,ysize=400 nsize = 500 ;IF (f1right*f1left GT 0) THEN BEGIN ; IF (f1right GT 0) THEN BEGIN ; f1left = -1.e-9 ; f1right = f1right ; ENDIF ELSE BEGIN ; f1left = f1left ; f1right = 1.e-9 ; ENDELSE ;ENDIF ;xrange = [f1left,f1right] ;f1 = dblarr(nsize) ;f1 = dindgen(nsize)/99.*(f1right-f1left) + f1left angle=dblarr(nsize) if n_elements(ks_range_deg) eq 0 then ks_range_deg=[-85,85] if ks_range_deg[1] gt 89.9 or $ ks_range_deg[1] gt 89.9 then begin message, 'ks_range_deg is outside of allowed range [-89.9,89.9]', /cont stop endif angle = dindgen(nsize)/float(nsize-1) * $ (ks_range_deg(1)-(ks_range_deg(0))) + ks_range_deg(0) angle = angle*3.14159/180.d0 f1 = tan(angle) ;translates to slopes from -11 to 11 xrange = [min([f1left,f1(0)]),max([f1right,f1(n_elements(f1)-1)])] ; Search kappa space, where kappa is only a function of f1 f2calc = dblarr(nsize) kappa = dblarr(nsize) high = 1.e+30 for ikl=0,nsize-1 do begin ; f1 loop f2calc(ikl) = f2func(f1(ikl),x,y,sigx,sigy,iswitch) kappa(ikl)=total((f1(ikl)*x+f2calc(ikl)-y)^2/$ (sigy^2+(f1(ikl)*sigx)^2) ) IF (kappa(ikl) lt high) then begin save1_kappa=ikl high = kappa(ikl) endif endfor print, 'kappa a minimum at (f1,f2):',f1(save1_kappa),f2calc(save1_kappa) print, 'kappa equals: ',kappa(save1_kappa) plot,f1,kappa,xtitle='f1',ytitle='kappa',xstyle=1,ystyle=1,xrange=xrange,title=titstr,/ylog oplot, [f1(save1_kappa),f1(save1_kappa)],[min(kappa),max(kappa)],line=0 time_stamp GOTO, zero_in ;=================BEGIN KAPPA SPACE SEARCH (NECESSARY SEARCH)================ search_kappaspace: IF ((verbose OR kappa_picture) AND NOT force_kappa_search) THEN print, 'kappa search necessary' ELSE IF ((verbose OR kappa_picture) AND force_kappa_search) THEN print, 'kappa search forced' IF (NOT print AND kappa_picture) THEN window,n_picture,xsize=550,ysize=400 nsize = 500 kappa_data = {f1:dblarr(nsize),kappa:dblarr(nsize),$ kappamin_index:0,$ both68:0,$ both95:0} angle=dblarr(nsize) angle = dindgen(nsize)/float(nsize-1)*$ (ks_range_deg(1)-ks_range_deg(0)) + ks_range_deg(0) angle = angle*3.14159/180.d0 f1 = tan(angle) ;translates to slopes from -11 to 11 xrange = [f1(0),f1(n_elements(f1)-1)] ; Search kappa space, where kappa is only a function of f1 f2calc = dblarr(nsize) kappa = dblarr(nsize) high = 1.e+30 save1_kappa = 0 for ikl=0,nsize-1 do begin ; f1 loop f2calc(ikl) = f2func(f1(ikl),x,y,sigx,sigy,iswitch) kappa(ikl)=total((f1(ikl)*x+f2calc(ikl)-y)^2/$ (sigy^2+(f1(ikl)*sigx)^2) ) if(kappa(ikl) lt high) then begin save1_kappa=ikl high = kappa(ikl) endif ENDFOR IF (verbose OR kappa_picture) THEN BEGIN print, 'kappa a minimum at (f1,f2):',f1(save1_kappa),f2calc(save1_kappa) print, 'kappa equals: ',kappa(save1_kappa) ENDIF IF (kappa_picture) THEN BEGIN plot,f1,kappa,xtitle='f1',ytitle='kappa',xstyle=1,ystyle=1,xrange=xrange,title=titstr,/ylog,yrange=[.1,100] oplot, [f1(save1_kappa),f1(save1_kappa)],[min(kappa),max(kappa)],line=0 time_stamp ENDIF kappa_data.f1 = f1 kappa_data.kappa = kappa kappa_data.kappamin_index = save1_kappa ; Need to bracket the zero knowing where this minimum is... ; First, do we have points on either side of the minimum? IF (save1_kappa EQ 0 OR save1_kappa EQ nsize-1) THEN BEGIN IF (verbose) then print, 'kappa search found no minimum...' f1new = slope_init f2new = intercept_init ireturncode = 2 GOTO, terminal ENDIF ; If we're here, the kappa search has found a minimum, and ; we're going to check and see what the derivatives look like f1dum = f1(save1_kappa) f2dum = f2calc(save1_kappa)*iswitch testmin = dkappadf1(f1dum,f2dum,x,y,sigx,sigy,iswitch) IF (testmin LT 0) THEN BEGIN f1dum = f1(save1_kappa+1) f2dum = f2calc(save1_kappa+1)*iswitch testdum = dkappadf1(f1dum,f2dum,x,y,sigx,sigy,iswitch) IF (testdum GT 0) THEN BEGIN IF (verbose OR kappa_picture) THEN print,'kappa search bracketed a minimum' f1left = f1(save1_kappa) f1right = f1(save1_kappa+1) testleft = testmin testright = testdum GOTO, zero_in ENDIF ELSE BEGIN IF(verbose OR kappa_picture) THEN print,'kappa search found a min and max between initial range' f1new = slope_init f2new = intercept_init ireturncode = 2 GOTO,terminal ENDELSE ENDIF ELSE BEGIN f1dum = f1(save1_kappa-1) f2dum = f2calc(save1_kappa-1)*iswitch testdum = dkappadf1(f1dum,f2dum,x,y,sigx,sigy,iswitch) IF (testdum LT 0) THEN BEGIN IF (verbose OR kappa_picture) THEN print,'kappa search bracketed a minimum' f1left = f1(save1_kappa-1) f1right = f1(save1_kappa) testleft = testdum testright = testmin GOTO, zero_in ENDIF ELSE BEGIN IF (verbose OR kappa_picture) THEN print,'kappa search found a min and max between initial range' f1new = slope_init f2new = intercept_init ireturncode = 2 GOTO,terminal ENDELSE ENDELSE ;==============END KAPPA SPACE SEARCH============================ ; Bolzano search for the zero of dkappadf1 proceeds zero_in: IF (verbose OR kappa_picture) THEN print, 'Zeroing in...' loop_counter = 0 loop: f1new = f1left- testleft*(f1right-f1left)/(testright-testleft) ;; jump_in: loop_counter = loop_counter+1 f2new=-total((f1new*x-y)/((f1new*sigx)^2+sigy^2))/$ total (1.0d0/((f1new*sigx)^2+sigy^2))$ *iswitch ;;ok jds f1dum=f1new f2dum=f2new*iswitch testnew = dkappadf1(f1dum,f2dum,x,y,sigx,sigy,iswitch) IF (verbose OR kappa_picture) THEN BEGIN print, loop_counter print, f1left,f1new,f1right,format='(3(e25.18,2x))' print, testleft,testnew,testright,format='(3(e25.18,2x))' print, '' ENDIF IF (kappa_picture) THEN BEGIN xyouts,f1new,.95*min(kappa),'+',alignment=.5 ; wait,.25 ENDIF ; Is testnew too close to 'zero' to pass up? IF (abs(testnew) LT zero) THEN BEGIN ireturncode = 0 goto, terminal ENDIF ; Have we been in this zero-in loop too long??? if (loop_counter ge n_iterations) then begin ireturncode = 5 goto, terminal endif ; Are f1left and f1right both close enough to f1new? IF (abs(f1new -f1left) LT (zero+abs(errorpercent/100.*f1left)) AND $ abs(f1new-f1right) LT (zero+abs(errorpercent/100.*f1right))) THEN BEGIN IF (testleft LT 0) THEN BEGIN ; Make sure it's a minimum ireturncode = 1 GOTO, terminal ENDIF ELSE BEGIN ireturncode = 3 GOTO, terminal ENDELSE ENDIF ; Is f1left close enough to f1new, but f1right isn't? IF (abs(f1new -f1left) LT (zero+abs(errorpercent/100.*f1left))) THEN BEGIN IF (verbose OR kappa_picture) THEN print, 'Flipping...' f1new=0.5*(f1left+f1right) GOTO, jump_in ENDIF ; Is f1right close enough to f1new, but f1left isn't? IF (abs(f1new-f1right) LT (zero+abs(errorpercent/100.*f1right))) THEN BEGIN IF (verbose OR kappa_picture) THEN print, 'Flipping...' f1new=0.5*(f1left+f1right) GOTO, jump_in ENDIF ; testnew and testright on same side of zero, re-bracket and try again IF (testnew*testright GT 0) THEN BEGIN f1right=f1new f2right=f2new testright = testnew ENDIF ; testnew and testleft on same side of zero, re-bracket and try again IF (testnew*testleft GT 0) THEN BEGIN f1left=f1new f2left=f2new testleft = testnew ENDIF GOTO, loop terminal: slope = f1new intercept = f2new ; will be zero for 1-param fit status = ireturncode IF (status EQ 1 OR status EQ 4) THEN BEGIN ; Evaluate chi2 and rchi2 s2= (slope*x+intercept-y)^2/(sigy^2+(slope*sigx)^2) chisq=total( s2 ) rchisq=chisq/double(npts-nparms-1) IF (rchisq LT 1 AND (verbose OR kappa_picture)) THEN print, 'rchi2 is < 1; RENORMALIZING chi2 surface (affects c68, c95, c99 and ks_c68, ks_c95, ks_c99)' ; Find errors on slope using kappa search method IF (ks_required AND status EQ 1) THEN BEGIN IF (rchisq LT 1) or keyword_set(renormalize) THEN BEGIN ; Must renormalize kappa surface kappa_temp = kappa/kappa(save1_kappa)*double(npts-2-1) ENDIF ELSE BEGIN kappa_temp = 0 ENDELSE find_ks_con_errors, ks_c68, ks_c95, ks_c99, nparms, f1, $ save1_kappa, kappa, kappa_data, nsize, kappa_picture, $ kappa_temp=kappa_temp, titstr=titstr ENDIF ; Find errors on slope and intercept using analytic method find_con_errors, c68, c95, c99, f1new, sigx, sigy, x, f2new, y, $ nparms, iswitch, status, rchisq IF (verbose OR kappa_picture) THEN BEGIN print, 'slope,intercept',slope,intercept print, 'rchisq',rchisq print, 'c68=[slope error,intercept error]=['+ $ string(c68(0),'(e9.2)')+','+string(c68(1),'(e9.2)')+']' print, 'c95=[slope error,intercept error]=['+ $ string(c95(0),'(e9.2)')+','+string(c95(1),'(e9.2)')+']' print, 'c99=[slope error,intercept error]=['+ $ string(c99(0),'(e9.2)')+','+string(c99(1),'(e9.2)')+']' IF (ks_required) THEN BEGIN print, 'ks_c68=[slope error left,slope error right]=['+ $ string(ks_c68(0),'(e9.2)')+','+string(ks_c68(1),'(e9.2)')+']' print, 'ks_c95=[slope error left,slope error right]=['+ $ string(ks_c95(0),'(e9.2)')+','+string(ks_c95(1),'(e9.2)')+']' print, 'ks_c99=[slope error left,slope error right]=['+ $ string(ks_c99(0),'(e9.2)')+','+string(ks_c99(1),'(e9.2)')+']' ENDIF print, 'status',status print, 'errorpercent',errorpercent ENDIF ENDIF ELSE BEGIN IF (verbose OR kappa_picture) THEN BEGIN print, 'Failed' print, 'status',status ENDIF ENDELSE IF (print) THEN BEGIN end_of_prog,/print ENDIF return END