00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 c------------------------------------------------------------------------
00021 c
00022 c Least-squares fitting routines.
00023 c
00024 c------------------------------------------------------------------------
00025
00026 subroutine lsfit1(x,y,n,w,nw,a,b)
00027 implicit real (a-h,o-z)
00028 save
00029 c
00030 dimension x(1),y(1),w(1)
00031 common/resid/resi
00032 c
00033 c Return least-squares best linear fit to y(x) (n points,
00034 c weighted by w) of the form y = a*x + b, where a is
00035 c specified and b is variable.
00036 c
00037 if(n.le.1)then
00038 b=0.
00039 resi=0.
00040 return
00041 end if
00042 c
00043 s1=0.
00044 s2=0.
00045 do i=1,n
00046 c
00047 if (nw.eq.1) then
00048 ww = w(1)
00049 else
00050 ww = w(i)
00051 end if
00052 c
00053 s1=s1+ww
00054 s2=s2+ww*(y(i)-a*x(i))
00055 end do
00056 b=s2/s1
00057 c
00058 resi=0.
00059 do i=1,n
00060 c
00061 if (nw.eq.1) then
00062 ww = w(1)
00063 else
00064 ww = w(i)
00065 end if
00066 c
00067 resi=resi+ww*(y(i)-a*x(i)-b)**2
00068 end do
00069 resi=sqrt(resi/s1)
00070 c
00071 dy=a*(x(1)-x(n))
00072 if(dy.ne.0.)then
00073 resi=resi/abs(dy)
00074 else
00075 resi=1.
00076 end if
00077 c
00078 end
00079
00080
00081 subroutine lsfit2(x,y,n,w,nw,a,b)
00082 implicit real (a-h,o-z)
00083 save
00084 c
00085 dimension x(1),y(1),w(1)
00086 common/resid/resi
00087 c
00088 c Return least-squares best linear fit to y(x) (n points,
00089 c weighted by w) of the form y = a*x + b, where a and b
00090 c are both variable.
00091 c
00092 if(n.le.1)then
00093 a=0.
00094 b=0.
00095 resi=0.
00096 return
00097 end if
00098 c
00099 c We will used NORMALIZED numbers for the fitting, to avoid
00100 c problems with single-precision roundoff...
00101 c
00102 xo = x(1)
00103 xfac = 1./(x(n) - xo)
00104 yo = y(1)
00105 yfac = 1./(y(n) - yo)
00106 c
00107 s1=0.
00108 sx=0.
00109 sy=0.
00110 sxx=0.
00111 sxy=0.
00112 do i=1,n
00113 c
00114 if (nw.eq.1) then
00115 ww = w(1)
00116 else
00117 ww = w(i)
00118 end if
00119 c
00120 xx = (x(i) - xo)*xfac
00121 yy = (y(i) - yo)*yfac
00122 c
00123 s1=s1+ww
00124 sx=sx+ww*xx
00125 sy=sy+ww*yy
00126 sxx=sxx+ww*xx**2
00127 sxy=sxy+ww*xx*yy
00128 end do
00129 d=sxx*s1-sx**2
00130 c
00131 a=0.
00132 b=0.
00133 if(d.ne.0.)then
00134 d=1./d
00135 a=d*(s1*sxy-sx*sy)
00136 b=d*(sxx*sy-sx*sxy)
00137 end if
00138 c
00139 c Now undo the normalization...
00140 c
00141 b = (b-a*xfac*xo)/yfac + yo
00142 a = a*xfac/yfac
00143 c
00144 resi=0.
00145 do i=1,n
00146 c
00147 if (nw.eq.1) then
00148 ww = w(1)
00149 else
00150 ww = w(i)
00151 end if
00152 c
00153 resi=resi+ww*(y(i)-a*x(i)-b)**2
00154 end do
00155 resi=sqrt(resi/s1)
00156 c
00157 dy=a*(x(1)-x(n))
00158 if(dy.ne.0.)then
00159 resi=resi/abs(dy)
00160 else
00161 resi=1.
00162 end if
00163 c
00164 end
