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 fitxy(input,istart,nin,x,y,z,w,n,iprompt,ier)
00021 save
00022 c
00023 c Perform least-squares fit to y(x).
00024 c NOTE: It is assumed that x is sorted!
00025 c
00026 character*(*) input
00027 dimension x(*),y(*),z(*),w(*)
00028 c ***** For now, w really is of dimension 1. *****
00029 c
00030 common/resid/resi
00031 c
00032 c Equal weights for all points.
00033 c
00034 nw = 1
00035 go to 10
00036 c
00037 entry fitxyz(input,istart,nin,x,y,z,w,n,iprompt,ier)
00038 c
00039 c Allow the weights to be specified in the z array...
00040 c
00041 nw = n
00042 c
00043 10 if (istart.gt.nin) then
00044 i1=1
00045 nn=n
00046 call minmax(x,n,x1,x2)
00047 else
00048 call readrq(input(istart:nin),2,
00049 & x1,x2,dum,dum,*1001)
00050 dx = x(2) - x(1)
00051 i1 = 0
00052 do 175 i=1,n
00053 c
00054 if (i.gt.1.and.(x(i)-x(i-1))*dx.lt.0.) then
00055 c
00056 c The x-array is not monotonic.
00057 c
00058 write(6,*)'The x array must be ordered.'
00059 & ' Use "so2 x y" to sort x and reorder y.'
00060 go to 1001
00061 end if
00062 c
00063 if ((x2-x(i))*(x(i)-x1).ge.0.) then
00064 if (i1.eq.0) i1 = i
00065 else if (i1.gt.0) then
00066 go to 176
00067 end if
00068 c
00069 175 continue
00070 i=n+1
00071 176 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 c
00078 x1 = x(i1)
00079 x2 = x(i1+nn-1)
00080 end if
00081 c
00082 if (nw.gt.1) then
00083 c
00084 c Make sure the weights are non-negative.
00085 c
00086 do 177 i=1,n
00087 if (z(i).lt.0.) nw = 1
00088 177 continue
00089 end if
00090 c
00091 if (nw.eq.1) then
00092 call lsfit2(x(i1),y(i1),nn,w(1),1,aa,bb)
00093 c ***** For now, w really is of dimension 1. *****
00094 else
00095 call lsfit2(x(i1),y(i1),nn,z(i1),nn,aa,bb)
00096 end if
00097 c
00098 c Print some basic data on the fit.
00099 c --------------------------------
00100 c
00101 if (iprompt.eq.1) then
00102 call devoff
00103 write(6,178)aa,bb,-bb/aa,resi
00104 178 format(' slope =',1pe11.3,', intercepts =',2e11.3,
00105 & ', residual =',e11.3)
00106 end if
00107 c
00108 if (input(3:3).eq.'p') then
00109 c
00110 c Plot the fit.
00111 c
00112 call fr inches(x1,aa*x1+bb,rr,ss)
00113 call devon
00114 call plotin(rr,ss,3)
00115 call fr inches(x2,aa*x2+bb,rr,ss)
00116 call plotin(rr,ss,2)
00117 end if
00118 c
00119 ier = 0
00120 return
00121 c
00122 1001 ier = 1
00123 return
00124 c
00125 end
