15.3 Straight-Line Data with Errors in Both Coordinates 669 sx=vector(1,ndat); sy=vector(1,ndat); ww-vector(1,ndat); avevar(x,ndat,&dumi,&varx); Find the x and y variances,and scale avevar(y,ndat,&dumi,&vary); the data into the global variables scale=sqrt(varx/vary); for communication with the func- nn=ndat; tion chixy. for (j=1;j<=ndat;j++){ xx[i]=x[i]; yy[j]=y[j]*scale; sx[j]=sigx[j]; sy[j]=sigy [j]*scale; ww[j]=sqrt(SQR(sx[j])+SQR(sy[j])); Use both z and y weights in first trial fit. fit(xx,yy,nn,ww,1,&dum1,b,&dum2,&dum3,&dum4,&dum5);Trial fit for b. 83g offs=ang[1]=0.0; Construct several angles for refer- 18881892 ang[2]=atan(*b): ence points,and make b anan- ang[4]=0.0; gle. 1800 ang [5]-ang[2]; ang[6]-POTN; for (j=4;i<=6;++)ch[i]=chixy(ang[i]); Cambridge from NUMERICAL RECIPESI mnbrak (&ang [1],&ang[2],&ang [3],&ch [1],&ch [2],&ch [3],chixy); Bracket the x-minimum and then locate it with brent. *chi2=brent (ang[1],ang [2],ang[3],chixy,ACC,b); *chi2=chixy(*b); server computer, (Nort 令 asaa; Compute x2 probability. America one paper University Press. THE *q=gammg(0.5*(nn-2),*ch12*0.5); for(r2=0.0,j=1;j<=nn;j++)r2+=ww[j]: Save the inverse sum of weights at ART r2=1.0/r2; the minimum bmx=BIG; Now,find standard errors for b as bmn=BIG; points where△x2=1. Programs offs=(*ch12)+1.0; for(j=1;j<=6;j++)[ Go through saved values to bracket if (ch[i]>offs){ the desired roots.Note period- d1=fabs(ang[j]-(*b)); icity in slope angles. while (d1 >PI)d1 -=PI; to dir d2=PI-d1: if (ang[j]*b){ swap=d1; d1=d2; ectcustser d2=swap; if (d1 bmx)bmx=d1; if (d2 bmn)bmn=d2; 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 10-621 -43108 if (bmx BIG){ Call zbrent to find the roots. bmx=zbrent (chixy,*b,*b+bmx,ACC)-(*b); amx=aa-(a); (outside bmn=zbrent (chixy,*b,*b-bmn,ACC)-(*b); amn=aa-(*a); North Software. *sigb=sgrt (0.5*(bmx*bmx+bmn*bmn))/(scale*SQR(cos(*b))); *siga=sqrt(0.5*(amx*amx+amn*amn)+r2)/scale; Error in a has additional piece Ame else (*sigb)=(*siga)=BIG; r2. a /scale: Unscale the answers. *b=tan(*b)/scale; free_vector(ww,1,ndat); free_vector(sy,1,ndat); free_vector(sx,1,ndat); free_vector(yy,1,ndat); free_vector(xx,1,ndat);15.3 Straight-Line Data with Errors in Both Coordinates 669 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). sx=vector(1,ndat); sy=vector(1,ndat); ww=vector(1,ndat); avevar(x,ndat,&dum1,&varx); Find the x and y variances, and scale the data into the global variables for communication with the func￾tion chixy. avevar(y,ndat,&dum1,&vary); scale=sqrt(varx/vary); nn=ndat; for (j=1;j<=ndat;j++) { xx[j]=x[j]; yy[j]=y[j]*scale; sx[j]=sigx[j]; sy[j]=sigy[j]*scale; ww[j]=sqrt(SQR(sx[j])+SQR(sy[j])); Use both x and y weights in first } trial fit. fit(xx,yy,nn,ww,1,&dum1,b,&dum2,&dum3,&dum4,&dum5); Trial fit for b. offs=ang[1]=0.0; Construct several angles for refer￾ence points, and make b an an￾gle. ang[2]=atan(*b); ang[4]=0.0; ang[5]=ang[2]; ang[6]=POTN; for (j=4;j<=6;j++) ch[j]=chixy(ang[j]); mnbrak(&ang[1],&ang[2],&ang[3],&ch[1],&ch[2],&ch[3],chixy); Bracket the χ2 minimum and then locate it with brent. *chi2=brent(ang[1],ang[2],ang[3],chixy,ACC,b); *chi2=chixy(*b); *a=aa; *q=gammq(0.5*(nn-2),*chi2*0.5); Compute χ2 probability. for (r2=0.0,j=1;j<=nn;j++) r2 += ww[j]; Save the inverse sum of weights at r2=1.0/r2; the minimum. bmx=BIG; Now, find standard errors for b as points where ∆χ bmn=BIG; 2 = 1. offs=(*chi2)+1.0; for (j=1;j<=6;j++) { Go through saved values to bracket the desired roots. Note period￾icity in slope angles. if (ch[j] > offs) { d1=fabs(ang[j]-(*b)); while (d1 >= PI) d1 -= PI; d2=PI-d1; if (ang[j] < *b) { swap=d1; d1=d2; d2=swap; } if (d1 < bmx) bmx=d1; if (d2 < bmn) bmn=d2; } } if (bmx < BIG) { Call zbrent to find the roots. bmx=zbrent(chixy,*b,*b+bmx,ACC)-(*b); amx=aa-(*a); bmn=zbrent(chixy,*b,*b-bmn,ACC)-(*b); amn=aa-(*a); *sigb=sqrt(0.5*(bmx*bmx+bmn*bmn))/(scale*SQR(cos(*b))); *siga=sqrt(0.5*(amx*amx+amn*amn)+r2)/scale; Error in a has additional piece } else (*sigb)=(*siga)=BIG; r2. *a /= scale; Unscale the answers. *b=tan(*b)/scale; free_vector(ww,1,ndat); free_vector(sy,1,ndat); free_vector(sx,1,ndat); free_vector(yy,1,ndat); free_vector(xx,1,ndat); }