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 subroutine genfit(input,istart,nin,x,y,z,n,nz,iprompt,*)
00021 save
00022 c
00023 c Perform a general least-squares fit to y(x).
00024 c
00025 character*(*) input
00026 c
00027 parameter (MMAX = 21)
00028 dimension x(1),y(1),z(1),a(MMAX)
00029 c
00030 external poly,trig
00031 c
00032 c Equal weights for all points if nz <= 0.
00033 c Otherwise, specify the weights in the z array...
00034 c
00035 nw = nz
00036 if (nz.le.0) nw = 1
00037 c
00038 c Determine the array ranges. Note that the arrays MUST be sorted
00039 c in x for this to work.
00040 c
00041 call sort3(n,x,y,z)
00042 c
00043 10 if (istart.gt.nin) then
00044 i1 = 1
00045 nn = n
00046 ma = 2
00047 else
00048 ma = 1
00049 x1 = -1.e20
00050 x2 = 1.e20
00051 call readrq(input(istart:nin),3,
00052 & xma,x1,x2,dum,*1001)
00053 ma = nint(xma)
00054 if (ma.le.0) ma = 1
00055 c
00056 ma = ma + 1
00057 if (ma.gt.MMAX) then
00058 if (iprompt.eq.1) then
00059 call devoff
00060 write(6,*)' Must specify order <= ',MMAX-1
00061 go to 1001
00062 end if
00063 end if
00064 c
00065 i1=0
00066 do 150 i=1,n
00067 if (x(i).gt.x1.and.i1.eq.0) i1=i
00068 if (x(i).gt.x2) go to 160
00069 150 continue
00070 i = n+1
00071 160 i1 = max(1,i1)
00072 nn = i-i1
00073 if (nn.le.1) then
00074 write(6,*)'No points in specified range.'
00075 go to 1001
00076 end if
00077 end if
00078 c
00079 c Note: Input ma is the order of the polynomial fit (<=20).
00080 c Internal ma is one greater than this number.
00081 c
00082 c Rescale the x and y data to lie in [-1,1], to reduce the effects of
00083 c finite-precision arithmetic. This appears to do a better job than
00084 c scaling to [0,1].
00085 c
00086 call minmax(x,n,xmin,xmax)
00087 xscale = 2./(xmax - xmin)
00088 xref = .5*(xmin+xmax)
00089 c
00090 call minmax(y,n,ymin,ymax)
00091 yscale = 2./(ymax - ymin)
00092 yref = .5*(ymin+ymax)
00093 c
00094 do 200 i=1,n
00095 x(i) = (x(i) - xref)*xscale
00096 y(i) = (y(i) - yref)*yscale
00097 200 continue
00098 c
00099 xscale = 1./xscale
00100 yscale = 1./yscale
00101 c
00102 x1 = x(i1)
00103 x2 = x(i1+nn-1)
00104 c
00105 c Use an intermediate C routine to allocate workspace space.
00106 c (Pointers are NOT standard FORTRAN yet!)
00107 c
00108 if (input(2:2).eq.'t') then
00109 call cfit(x(i1),y(i1),z(i1),nn,a,ma,trig,nw,chisq,iret)
00110 else
00111 call cfit(x(i1),y(i1),z(i1),nn,a,ma,poly,nw,chisq,iret)
00112 end if
00113 c
00114 c Print out some basic data on the fit obtained.
00115 c
00116 if (iprompt.eq.1) then
00117 call devoff
00118 write(6,250)(a(i),i=1,ma)
00119 250 format(' Coefficients: ',1p,5e11.3,:/(15x,5e11.3))
00120 write(6,*)' Chisq = ',chisq
00121 end if
00122 c
00123 if (input(3:3).eq.'p') then
00124 c
00125 c Plot the fit.
00126 c
00127 call uplotin(x1*xscale+xref,
00128 & afunc(x1,a,ma,input(2:2))*yscale+yref,3)
00129 c
00130 do 300 i=1,n
00131 300 if (x(i).gt.x1.and.x(i).le.x2)
00132 & call uplotin(x(i)*xscale+xref,
00133 & afunc(x(i),a,ma,input(2:2))*yscale+yref,2)
00134 end if
00135 c
00136 c Unscale x and y before returning.
00137 c
00138 do 350 i=1,n
00139 x(i) = xscale*x(i) + xref
00140 y(i) = yscale*y(i) + yref
00141 350 continue
00142 c
00143 return
00144 1001 return 1
00145 c
00146 end
00147
00148
00149 function afunc(x,a,ma,which)
00150 save
00151 c
00152 parameter (MMAX = 21)
00153 dimension a(ma),f(MMAX)
00154 character*1 which
00155 c
00156 if (which.eq.'t') then
00157 call trig(x,f,ma)
00158 else
00159 call poly(x,f,ma)
00160 end if
00161 c
00162 sum = 0.
00163 do 10 i=1,ma
00164 10 sum = sum + a(i)*f(i)
00165 c
00166 afunc = sum
00167 c
00168 return
00169 end
