subroutine lsqfit(x,y,npts,axlen,iflag) C C subroutine to generate a least squares fitted C line through the co-ordinate data pair C C x,y(real)=data arrays C npts(int)=number of points in the array C note:array(npts+1)should contain the minimum C array(npts+2)should contain the increment C axlen(real)=length in inches of the x and y axes. C The parameters will be plotted in three lines,using the C first inch above the y-axis. C iflag(int)=flag to suppress/plot statistical parameters: C 0=suppress plotting of parameters C 1=plot values for slope,intercept and correlation C real x(1),y(1) double precision xvar,yvar,sumx,sumy,sumx2,sumy2,prodxy,xmean,ymean double precision interc,interr,slope,sloper,correl,ydev,xdev,smean double precision dnpts idec=8 !idec=the no. of decimal places to print C C Note: Rather than using the DATA statement to initialize variables, C assignment statements have been used.There appears to be a C problem with the DATA statements when the subroutine is called C within a loop,such that the variables are initialized only with C the first pass of the loop, and subsequent iterations will give C erroneous results. C sumx=0D0 sumy=0D0 prodxy=0D0 sumx2=0D0 sumy2=0D0 smean=0D0 dnpts=dble(float(npts)) do 10 i=1,npts sumx=sumx+x(i) !find sum of x array sumy=sumy+y(i) !find sum of y array prodxy=prodxy+x(i)*y(i) !find sum of product of x*y sumx2=sumx2+x(i)**2 !find sum of x squared sumy2=sumy2+y(i)**2 !find sum of y squared 10 continue xmean=sumx/dnpts !find mean of x ymean=sumy/dnpts !find mean of y slope=(xmean*sumy-prodxy)/(xmean*sumx-sumx2) interc=ymean-slope*xmean do 20 i=1,npts !variance in y 20 smean=smean+(x(i)-xmean)**2 yvar=sumy2/dnpts-ymean**2 xvar=sumx2/dnpts-xmean**2 ydev=dsqrt(yvar) !standard deviation in y xdev=dsqrt(xvar) correl=slope*xdev/ydev sloper=ydev/sqrt(smean) interr=ydev*sqrt(1d0/dnpts+xmean**2/smean) C C Now calculate the end points of the fitted line C xmin=x(npts+1) dx=x(npts+2) ymin=y(npts+1) dy=y(npts+2) ymax=ymin+dy*axlen xmax=xmin+dx*axlen ypos=xmin*slope+interc xpos=xmin xpt=xmax ypt=xmax*slope+interc C C clipping routine C if(ypt.gt.ymax)xpt=(ymax-interc)/slope if(ypt.lt.ymin)xpt=(ymin-interc)/slope if(ypt.gt.ymax)ypt=ymax if(ypt.lt.ymin)ypt=ymin if(ypos.gt.ymax)xpos=(ymax-interc)/slope if(ypos.lt.ymin)xpos=(ymin-interc)/slope if(ypos.gt.ymax)ypos=ymax if(ypos.lt.ymin)ypos=ymin C C scale the points C xpos=(xpos-xmin)/dx xpt=(xpt-xmin)/dx ypos=(ypos-ymin)/dy ypt=(ypt-ymin)/dy C C now plot the line C call plot(xpos,ypos,3) !move to the first point call plot(xpt,ypt,2) !draw the line to the last point if(iflag.eq.0)return call symbol(.25,axlen+.25,.14,12HCorrelation:,0.,12) call number(2.00,axlen+.25,.14,correl,0.,idec) call symbol(3.05,axlen+.25,.14,8hsigma(Y),0.,8) call number(4.15,axlen+.25,.14,ydev,0.,idec) call symbol(.25,axlen+.5,.14,10HIntercept:,0.,10) call number(2.00,axlen+.5,.14,interc,0.,idec) call symbol(3.15,axlen+.5,.14,3h+/-,0.,3) call number(4.15,axlen+.5,.14,interr,0.,idec) call symbol(.25,axlen+.75,.14,6HSlope:,0.,6) call number(2.00,axlen+.75,.14,slope,0.,idec) call symbol(3.15,axlen+.75,.14,3H+/- ,0.,3) call number(4.15,axlen+.75,.14,sloper,0.,idec) return end