LSQ DOCLSQF COMLSQF FOR1LSQW BAS?NBS3 DATRNDTEST5DAT À y(i) yc(i)', x ' dy(i)',a1) 180 format(1x,i3,1p4e15.7,a1) 190 formatP ,=LSQF1/R }+?Command error?File not found?Can't enter file~/#+!+,!+,!++Á,COMSCN,"U&,e,:4,G:3,O*.[,# M,e, [,D,2,*+++"%*U&:%2%*%t/:r+ ʁ,2%2%2%2%2%2%2%2%<2%2%,ů2% .a&-V& ,>2%!%-> 2V&2^&2U&y,2%,:%2% .g&:%-d&-V&:U&2% 2%$-:%@2%!%-y=7 . 7-U&L- 7!%-:%!%-!%:%|-%d&> =r--:%‘-:%‘-:%2% _2͊0+0A0*%##:%2%GU& ~--w# - #- w#-:%!^&~ .w#-S/:7>..x ./Ox=yO.[A@2U&S/x <. V&ʘ./C.Th./T7/Y7..L‚./S7/T7>..R7/D7R7>..> ک.6 #=.y.2%S/^&x../.y/y/L.2%/O.N.2%.O.2%.R/2%.C/2%2%.M!/2%.I,/2%.X9/>2%.ZE/2%.P7:%<2%./[q/0q/Am/:q/V/O!~6Gʣ/#~ Œ/~/q+x+#ʝ/~Ò/2%/2%:%¶/!0>*/ p+> /:q+!r+"+o>g6 #> w~# / / 0#"+_*++~ #t/#"++~7ȷ> /> /n'!"l'!%:%T0l&!"j&!%:%ʅ0~0ͧ0< ,ͧ0<,!"n)p)!%~0y0 w#°0w%0p):%-:%<0j1> |1> |1>|1*l'n'%T1n'%0:%-:%0>1*j&l&%T1l&%0|}_1ƀ)`ió1> |1> |1> |1:%-:%§2 =³1!l&>!%ó1:%-:%/*j&|=1#"j&k&w:%<02*n)|+2#"n)o)0p)>%=T2!=02!>P2_!% y2 w2 |1~€2Ñ27ȷ!2+@@¸2> /!2+!0DISK FULL!2:>’3@3w#2%.ʽ3!ڽ3w#–32%.3!ҫ336 #½3:%.3w#3p2232222:2=!42264!2"2>22*2~14764#"2ɆADD ADIANA ģANA ŦANA ANICALL/CM?CMøCMP ¾CMP CPIDAD )DAD +DCX =DCR MVI A.MVI LORA ȵORA ̶ORA PCHPOP PUSH RARARERLßSBB SBB ȞSBISUB SUI"SHLD2STASUB XCHǯXRA XRA īXRA ŮXRA XRI>oͫG>ͫG>ͫG>gͫGͥU*aT͈>:$7*aT2>9*aT=©6>>ͫG:G>ͫG>>ͫG:G>ͫG6>ͫG͈>:(A6>ͫG>{ͫG>ͫG>zͫG>ͫG>*"]>2e=2!b?**udDM!"*"x7w# x7"!"N*aT:=*aT>2Z7>9>ͫG*aT:6:=ʆ7>9>GͫG*aT>2͈>>ͫG>xͫG6:=a7*aTW?>ͫG*aT>2=9>6!"7ͳA͑!:7*" ""ß%@2 :7ͳAoG>H8{%xͳA>ͫG>2s:,+8:ʟ% ͳA2sgo9"S>(-h͑!:sK8f*$8:,38v\:s8&hn%͈>:(s8A62(GC͜:%8c]kͷO!8nG*+V+^+aT:%»8c]kͷO!8n/P̓?=g˜8* 2(1X,'Too Many Parameters !!!')9*: 9 (">/ͫG><ͫGW9*a(1X,'Least Squares Polynomial'/' Double Precision'/' (c) 9/ 8/1986 R. H. Sherman'/' 50 points maximum'/)(' Is Data input from Keyborad (K) or File (F) ?')(11A1)(1X,'Input Degree of Polynomial to be fit = ')(I1)(' Data filename'/' 12345678901'/)('+Input drive ID - 0, 1, 2 :')(10A8)(I2)(2F12.0)(/1X,10A8)(1X,I3,' data points read'/)('+Input number of data points = ')(1X,'Enter Data:'/)('+ x(',I2,') = ')('+ y(',I2,') = ')(1X,'Data:')(1X,I3,1X,1P2G13.6)(/1X,'Matrix Coefficients:')(1X,I3,1X,1P6G14.6/18X,6G14.6)(/1X,'Inverse Matrix Coefficients:')(/1X,'Fitted Coefficients and std. deviations:')(1X,I3,2X,1PG20.10,2X,E15.7)(/' i x(i) y(i) yc(i)',' dy(i)')(1X,I3,1P4E15.7)(/1X,'Wvar, Standard Deviation of the fit = ',1P2E16.7)(1X,'Create an output file - Y/N ? ')(' Run another case - Y/N ?')B2IB!1,4!1,4!.,G!>ͻ+4*G}2!K "G*G}2!1,4!.,K!>+4!K͜ 7!1,4!.,!>ͻ+4K!ͽ?W!.,!>+4!"a!,*)))"K"*!>+*"*!>+4*#>3>_*+"Ei!1,!>+4s!1,E!>ͻ+4å!1,4]!.,E!>ͻ+4!1,4!"!1,!>ͻ+4a!.,*)))K"*!>+4!1,!>ͻ+4a!.,*)))"*!>+4*#*E{z!1,4]!.,C!>ͻ+4*C#"C*E*C{ozg}?2!1,4å!"*" *)))"* !Pt+* "!ͳ:ͩ)*&!""*)))"* !Pt+* "*")))K"* ++*"$$*!*#(*&*"#*E{z8* )))"*!Pt+* "*)))"* !Pt+* )*&* #*C{z*#*C{z!"*)))"!ͳ:ͩ)*&!" *)))"* )))""*K"$*+"&&*$!*+'*#(*&* #*E{z"*#*C{z!1,4!1,!"!>ͻ+*)))"K"*!>+*"*!>+*#*E{z4!1,4!"$!1,!>ͻ+*)))"*!>+!" *)))"* !Pt+* "*!>+* #*C{z;4*#*C{z*a!CWB!1,4!"$!1,!>ͻ+!" *)))"* !Pt+*) "*!>+* #*C{z4*#*C{z!"*)))"!ͳ:ͩ)*&!" *)))""* !Pt+*) "* ))))*+'*#(*&* #*C{z7*#*C{z!ͳ:ͩ)!.&!"*)))""*K"C* *&*k"$*)*(U**$&*)))k")*+'!.#(!.&*#*E{z*E*C{ozg}=2e!ͳ:ͩ)!3&Ç*E*C{ozg"!.)*=&!3&!3%!;&!"*)))""*!Pt+*) )!3+'!6&!6%*&*#*C{z*G}2i!1,!>+4f!1,4!"!1,!>ͻ+*)))""*!>+*"*!>+4*#*C{z$!1,4!"!1,!>ͻ+*)))"K"*!>+*"*!>+*"$*$!>+*k"&*&!>+4*#*E{z!1,;3!>+4=!1,4!.,>!>ͻ+4*>}2q!Y ">*>}2ʋ͏#b!1,4!.,@!>ͻ+4*@}2!Y "@*@}2¥!͢3͵*  b>22W>e?*a|aT:%F:*"](1X,'***** SINGULAR MATRIX *****')"">! t;!"!"*)))"*!Pt+** "! ͳ:ͩ)*&**{ozg}2*)))"*!Pt+** "! ͳ:ͩ)*&*#*~#fo{zo*#*~#fo{zi!"*"*)))"*!Pt+**)2‚*#*~#fo{z&5! 1,4!* s#r!"*)))"*!Pt+"*"!*"#)!&*)))"%**%*)*#&**%*"'!)*'&*!* )!&*!* ")*)))"+*!Pt+*+* )*)&*)))"*!Pt+** "!)*&*#*~#fo{z*)))"*!Pt+**"!͠)*K&!-&!"*)))"*!Pt+**")!-+'*&*)))"*!Pt+** ")!-+'*&*#*~#fo{z!"**{ozg}2T *)))"*!Pt+**)U*!-&!"*)))"*!Pt+"**"!*)))"#**#*)!-+'*!#(*!&*)))"*!Pt+"** "!*)))"#**#* )!-+'*!#(*!&*#*~#fo{z*#*~#fo{zA*#*~#fo{z !* s#rger Quant" !" * +* " * ͩ%* w* #> > referencIllegal State" " `i" * ~#fo+}2 «!* ~#fo)))* )* +'! &* ~#fo++}2 ˜!* ~#fo++" !" * ~#fo* {ozg" * )))* )! #(* +'! &* #* {zM!* )! #(! &* )! &oFunction Call "!"!!/"ͳ:ͩ)!!&*!~#fo}2!!!!)!"!*!)!!+'!!&*!#*!~#fo{z!!!)unDuplicate Statement Labe(' Data filename'/' 12345678901'/)(11A1)('+Input drive ID - 0, 1, 2 :')(I1)(/1X,10A8,A1)(/1X,'Fitted Coefficients and std. deviations:',A1)(1X,I3,2X,1PG20.10,2X,E15.7,A1)(/' i x(i) y(i) yc(i)',' dy(i)',A1)(1X,I3,1P4E15.7,A1)(/1X,'Wvar, Standard Deviation of the fit = ',1P2E16.7,A1)h> 2?":I=2@"#O"!%1,4q"!%.,4"!%>+4!4"͜ w"!%1,4"!%.,A"!%>ͻ+4A"4"!%ͽ?*G}2@"3$"!%1,!%>+?"!%>+4"!%1,?"!%>+4!"C""!%1,C"!%>ͻ+*C")))"E""G"*G"!%>+*E""I"*I"!%>+?"!%>+4*C"#*C{zM$"!%1,?"!%>+4!"C"A#!%1,C"!%>ͻ+*C")))"E"K"G"*G"!%>+*E""I"*I"!%>+*E""K"*K"!%>+*E"k"M"*M"!%>+?"!%>+4*C"#*E{z$T#!%1,;3!%>+?"!%>+4>2I sOG:%"%*%~2%*%~z=Ɓ2%*%~a?!%/2%%:%og!}2%:%)+&!C~>@w&>&&H&&0)!'+'&'=%!Cw&SA&ͳ:ͩ)E&&͠)͵&!C)ʱ&:C::C/F+N+=&!C'!C>q#p#= '!C+'"C&&&*C)W::CW:!CN#F*U)qh'Gܚ(͘)x Q'H'(!Ct)`'j0 TeB׳]h!I.k &ͳ:ͩ)ï'&͠)!C+'SA!Cq#~++w+q':DNn"~`35zr1{r1h!I):yO2C)()::Cʵ&;(/<͵&9::CO͏+G2C!Cl):C2Cx!CC(ͣ(r(4͘):C(!C~++w:!C4#‹(4ʵ'+6!CC# £(ɯ# °(Ͱ((/!COyw#(G:C(!CVwz# (x(W:!C)(xr(!CwW:r(~w# )&ͳ:ͩ)0)&͠)!C#(&ͳ:ͩ)L)&͠)͵&!C(y2C!CC~q+b)qڃ)NsY+x)n) W~w+)Æ)!CÊ)/+é)ͳ:!Cw#¯):)}))))͏+!CF#^#V#N7:)}C~#)G++Ny:C)!Cͳ:>8**SA:F*:4*!C6͜+g>P;|]+!> *akͳ:U**U*:!C~+>w:͓*> 2C!PA"C!US"C!E "Cͤ**B6 B*B6 T@BC~#š*!C ~B#©*͓*> 2C!ST"C!OP"C! "Cͤ*@B!C~Gx*ƀV:w͏+w+ɷ&:W::x{ҿ:!C6͜+wg+z+ >P;|I:{>2Cy9/))Ҋ+ =‚+!C~7w?##wy7O*C*C"C"Ca!+!+!+!+!"hC͆;͚;~#fo"PC͚;"&C*PC"*C:C,͘-+5+>2rC2sC<2TCG,>2rC6,2rC<6,>2rC2sC2TC~#fo"C!C"-C~#+U,!"VC!9##"+:sCy,go"XC"ZC"^C"\Cõ,!XC>t;*XC~#,go"XC*ZC~#,go"ZC*\C~#fo"\C*^C~#foµ,"^C2tC2C2C2C2cC2HCgo"`C"*Ca:*VC|,^6*`C"uC(ʴ5SA21C22C2/C20C<2C:rC>ï39]0- -> 6w -:JC/<0-0-6=)-2JC:hC7:iC/ 76=A-9:JCf-f- ͚4=Z-2JC:hCx->6=q-:jCG͚46x-!jCCw# Ž-:rC-*1CDM*/C"C*-C *&C-`i"1C"&C*/CDM*/D"C*&C*-C -`i"/Cü-:iC*C+}|&.:rC.*C"/C>ͯ3*-C->ͯ3*/C"C*-C-w#=-**C+"*C|-<G.>G!3C60#9!JC4:qC0w:M."fC~:ځ. w+q.w!3C|–.}–.!JC4#"MC|/W}/_*fC"fC>2HC:hCڝ1&8:ICG:C.x2GC?.:C.2LCL4:C.0͚4!GC5 /<.:GC / /'4.E:HC/D͚4:JC!C2JC)/+0/-/<2JC͚4:JC_! ;{0G͚4}0G͚4:hCڝ1&8:ICG!JC:Cwq/?.2LC:IC2GC:JC/2LC!GCwL4:JC//*fC"dC!"fC'4!JC4¤/*dC"fC'4ô/:hCڝ1&8:JC.!IC/<.Fw!RC~w2bC!C~w<2cC]/2C2RC2IC2cC>2bC ͚40*`C*VCDM*C͎0+|0"C"`C 7*`C*VCF͚4*`C#"`C*C+"C|#@0 7*1CDM!bC5ʄ0*/C+}|ډ0*-C ~2C`i"1C>z0> z0]0ʗ0 !bC4Î0:UCI0:hC=ʹ0==ʹ0=ʹ00:IC0SA!nC0SA0#0:C7:hC0j9K97>2,C:hC1SA!IC:KC1:SC1:OCw~/2qC&8:JC?.:JC2GC2LCL4!kC~-2/W+~/_>2C2JC! ;}0o!JC4{1!3C:JCG2GC2LC2fCs#2!3C"MCL4V+^1L3C2͎0 C2>00T2:T202C>ɯ2IC2KC2SCL3-r2>2Cw2+z2L37*1C+"1C!bC43:C.¢2>2KC2JC3:C0Eʲ2D·2L3+2-2>2,CL3y2SC*1C+"1C!bC4322Ox 2SA> G2:,Cx3/<2JCG!IC~wÞ092:qC33 : :ͅ- ::C!Cw#w83ͅ-!IC:JCw3͎0 L3:RC<2bCL31Fo3T1>2jC͎0s37!jC~ʇ3>T>Fw> #w#w#w>2hC2iCK-> é3> é3>O}2CO!::CMAMA=_^#V ^#V33*ZC|3SA*XC|3SA*+!CjCw#4!jCC3:C!,Cwʅ5:rC>ï3!GC584:dC*4!fC504*MCF#"MCÚ4!RC~!C!GC2dCj442GC/4:C-Ě4!LC54.͚4'4s4SA ͚4!dC5‹4j4*/C*1C{zڻ4#"1C*-C+p!bC5ÿ4SA:cC:C 75!tC~44SA *fC"uC*C"{C 7!"bC*`C*VC:rC#5~'5͎0w#5#~'5"`C 7~':5G͚4*`C#"`C##5*`C#"`C#~')5 7:tC=a5!C4!"`C52tC 7*{C+"{C| 7*uC"`C!tC4 7!,C6:rC©5:TC©5*1Cz¢5*-C6 "/C8:,C21C22C:rC5*-D"/CP6!"1C 7!C~5w5**C+"*C|5:C€5*C| 7+"C:OC2IC2C2JC2,Cgo"jC"lC"nC"pC:RC!bCw:rC-64*}C:iCG:hC*&CjCB6nC~#B6"&C)6:TCʯ3*-D"/Cɯ2,Cgo"jC*`C*VC~# o6-†6,C o6"`C2CG06 Ҩ6*jCT])))_d6*jC|³6#x:,C6o>g"jC|x!qC~w+ 66!qC~w+6!C66kC!jCw#76 ^6"C)7, 7/ʀ5)M5'4P77:jC2C 7|G7SA!"C2UCG!7b7# y xQ7SAyO!7:rCt7!7 ^#V"}Cyґ7*`C"fC!"bC^6›7SA}2RC2OCy ڻ7:C.7^6*jC}2OC*`C+"`C5SAXH(AILEFGD804-X2U3X2X2X2X28804K-1|3.R//.:iCG*&CjC:hC8nCw#8"&C*C| 75!qC~ʾ8x2,Cw:hC=T8==T8=W8+>2C>wf8j9:mC8K9>2C7:͖9!,C~wg889j8£8!CjC39j8!jCC3!JC4/G6!jC~#±8f8f82,C!jCw#8!JC4>ï3*C*1C"1C*-DDMyx 7SA */C"1C 7 ͚4*C+"C|8 7:iCG:RC2JC!,C>w!JC5636Å-!jC~/G#~/G97; 6;7>{_zW}o|g=-;|g}o_;CZQR; o-yOzW{_xGb;`iN#F#q#p#=w;2"C"wC"yC`i"#C2%C:%Cү;!wC_;:"C*#C;_~#fo:%C<2%CW:"C;h: 3fͥU*a*: eDMæe:eͧͳA:f*g|f"]1DFORTDAT>>i>@>>MA >SASASASA!;:C_~_<:C_!;^#V!;:C_~!<:C_^#Vz #w##<ͧ<>wnB%>|=}="<">B͖F< >6b?z::(>0Ͱ<*-C*/C&\>w#S>ͧ?7Ͱ<*-C>2/C͝?? Ͱ<2<:/C:/C>SA> ?*-C~H?#=¦>> H?F<w:<7Ͱ<2/C2< ??>?*/C&*-Cw"/C >}>SA:<?ɷF<?n<6n<4?5!;/Ɓ_#>WR<͝?2<>w7nWɇk?>H?n<6_<~2Ca:F ͗@ͬ@2?:?*?> ͗@> ×@*>B"@|E@}H@"?">BMA:?̤@:/C*-C=> ͗@~+ʉ@1y@> ͗@É@> ͗@~0‰@> ͗@#~#͗@=Ë@_>2?5@@@@@@@@MA2/CB @*/C&*-Cw"/C }@SA:/C*-C=> B~+&A1A> B&A> B~ &A0&A> B#~#B=(A **** at address ** SA!4A͝A~"~eA>#O!A ~B#~B:B<2BҊAҚA!9A͝AͧA͚A@B!HA~B#ÝA|ʹA}ʹA> BͽAƐ':'BIDF0MPIRFWITEXDOMLDZLGSQIBTLOBDEISBEINOVCNGLGSSNA2IODTBIRCEFFNDFUNOM??DFUNO2B2B!=B">B*+ 2B!:B!!""BUB!HB*>B> B> BP BP :B_*B~ B#~"XB_*B~ C#~"B VB!WB^#6 #> BwBw͖'mnnn+n܄ program lspoly c Double Precision Least squares Polynomial c with input and output files c 9/08/1986 c (c) Robert H. Sherman c 1931 Camino Manzana, Los Alamos, New Mexico 87544 c byte fnmi(11), fnmo(11), iopn implicit double precision(a-h,o-z) dimension x(50),y(50),dy(50),yc(50),xtitle(10) dimension a(10,10),b(10,10),c(10),r(10),sdc(10) common xtitle,x,y,dy,yc,c,sdc,wvar,sd,n,np,ipath,iopn c iopn=0 write(1,300) write(1,310) read(1,320)ipath if (ipath .eq. 1hk) ipath=1hK if (ipath .eq. 1hK ) go to 20 write (1,350) read (1,320) fnmi call upper (fnmi) write (1,360) read (1,340) id call open (10,fnmi,id) read(10,365)xtitle do 10 i=1,51 10 read (10,380,end=15) x(i),y(i) 15 np=i-1 write(1,381)xtitle write(1,382)np go to 40 20 write(1,390) read(1,370)np write(1,400) do 30 i=1,np write(1,410)i read(1,380)x(i) write(1,420)i 30 read(1,380)y(i) 40 write(1,330) read(1,370)n n=n+1 if (np .ge. n ) go to 45 write(1,421) 421 format(1x,'Too Many Parameters !!!') go to 40 45 do 60 i=1,n do 60 j=i,n a(i,j)=0. do 50 k=1,np 50 a(i,j)=a(i,j)+pwr(x(k),i+j-2) 60 a(j,i)=a(i,j) do 70 i=1,n r(i)=0. do 70 j=1,np 70 r(i)=r(i)+y(j)*pwr(x(j),i-1) write(1,430) write(1,440)(i,x(i),y(i),i=1,np) write(1,450) do 80 i=1,n 80 write(1,460)i,r(i),(a(i,j),j=1,n) c call matinv(n,a,b,ierr) c write(1,470) do 90 i=1,n 90 write(1,460)i,(b(i,j),j=1,n) do 100 i=1,n c(i)=0. do 100 j=1,n 100 c(i)=c(i)+b(i,j)*r(j) sumd=0. do 110 i=1,np yc(i)=poly(x(i),n,c) dy(i)=yc(i)-y(i) 110 sumd=sumd+dy(i)*dy(i) if (np .gt. n) go to 112 wvar=0. go to 113 112 wvar=sumd/(np-n) 113 sd=dsqrt(wvar) do 120 i=1,n 120 sdc(i)=dsqrt(wvar*b(i,i)) if (ipath .ne. 1hK ) write(1,381)xtitle write(1,480) do 130 i=1,n 130 write(1,490)i,c(i),sdc(i) write(1,500) do 140 i=1,np 140 write(1,510)i,x(i),y(i),yc(i),dy(i) write(1,520)wvar,sd c write(1,530) read(1,320)iprint if (iprint .eq. 1hy ) iprint=1hY if (iprint .eq. 1hY ) call output c write(1,540) read(1,320)iagain if (iagain .eq. 1hy) iagain=1hY if (iagain .eq. 1hY ) go to 40 endfile 6 stop c 300 format(1x,'Least Squares Polynomial'/' Double Precision'/ 1 ' (c) 9/ 8/1986 R. H. Sherman'/' 50 points maximum'/) 310 format(' Is Data input from Keyborad (K) or File (F) ?') 320 format (11a1) 330 format(1x,'Input Degree of Polynomial to be fit = ') 340 format (i1) 350 format (' Data filename'/' 12345678901'/) 360 format ('+Input drive ID - 0, 1, 2 :') 365 format(10a8) 370 format (i2) 380 format (2f12.0) 381 format(/1x,10a8) 382 format(1x,i3,' data points read'/) 390 format('+Input number of data points = ') 400 format(1x,'Enter Data:'/) 410 format('+ x(',i2,') = ') 420 format('+ y(',i2,') = ') 430 format(1x,'Data:') 440 format(1x,i3,1x,1p2g13.6) 450 format(/1x,'Matrix Coefficients:') 460 format(1x,i3,1x,1p6g14.6/18x,6g14.6) 470 format(/1x,'Inverse Matrix Coefficients:') 480 format(/1x,'Fitted Coefficients and std. deviations:') 490 format(1x,i3,2x,1pg20.10,2x,e15.7) 500 format(/' i x(i) y(i) yc(i)', x ' dy(i)') 510 format(1x,i3,1p4e15.7) 520 format(/1x,'Wvar, Standard Deviation of the fit = ',1p2e16.7) 530 format(1x,'Create an output file - Y/N ? ') 540 format(' Run another case - Y/N ?') c end subroutine matinv(n,a,b,ierr) c c Invert matrix c implicit double precision(a-h,o-z) dimension a(10,10),b(10,10) c do 10 i=1,n do 10 j=1,n b(i,j)=0. if ( i .ne. j ) go to 10 b(i,j)=1. 10 continue do 70 j=1,n do 20 i=j,n 20 if ( a(i,j) .ne. 0. ) go to 30 write(1,200) ierr=-1 return 30 do 40 k=1,n s=a(j,k) a(j,k)=a(i,k) a(i,k)=s s=b(j,k) b(j,k)=b(k,j) 40 b(i,k)=s t=1/a(j,j) do 50 k=1,n a(j,k)=t*a(j,k) 50 b(j,k)=t*b(j,k) do 70 i=1,n if ( i .eq. j ) go to 70 t=-a(i,j) do 60 k=1,n a(i,k)=a(i,k)+t*a(j,k) 60 b(i,k)=b(i,k)+t*b(j,k) 70 continue ierr=1 return c 200 format(1x,'***** SINGULAR MATRIX *****') c end subroutine upper (x) c c comvert filename to upper case c byte x(11), cupper do 10 i=1,11 10 x(i)=cupper(x(i)) return end function poly(x,n,c) c c evaluate polynomial (double precision) c implicit double precision(a-h,o-z) dimension c(10) if (n .eq. 1) go to 30 poly=c(n)*x if ( n .eq. 2 ) go to 20 nm2=n-2 do 10 i=1,nm2 j=n-i 10 poly=(poly+c(j))*x 20 poly=poly+c(1) return 30 poly=c(1) return end function pwr(x,n) c c x**n by multiplication (for double precision) c implicit double precision(a-h,o-z) pwr=1. if (n .eq. 0 ) return do 10 i=1,n 10 pwr=pwr*x return end subroutine output c byte fnmo(11),lf,iopn implicit double precision(a-h,o-z) dimension x(50),y(50),dy(50),yc(50),xtitle(10),c(10),sdc(10) common xtitle,x,y,dy,yc,c,sdc,wvar,sd,n,np,ipath,iopn c lf=x'0A' if (iopn .eq. 1 ) go to 10 write (1,100) read (1,110) fnmo call upper (fnmo) write (1,120) read (1,130) id call open (6,fnmo,id) 10 if (ipath.eq. 1hK ) go to 20 write(6,140)xtitle,lf 20 write(6,150)lf do 30 i=1,n 30 write(6,160)i,c(i),sdc(i),lf write(6,170)lf do 40 i=1,np 40 write(6,180)i,x(i),y(i),yc(i),dy(i),lf write(6,190)wvar,sd,lf iopn=1 return c 100 format (' Data filename'/' 12345678901'/) 110 format (11a1) 120 format ('+Input drive ID - 0, 1, 2 :') 130 format (i1) 140 format(/1x,10a8,a1) 150 format(/1x,'Fitted Coefficients and std. deviations:',a1) 160 format(1x,i3,2x,1pg20.10,2x,e15.7,a1) 170 format(/' i x(i) y(i) yc(i)', x ' dy(i)',a1) 180 format(1x,i3,1p4e15.7,a1) 190 format(/1x,'Wvar, Standard Deviation of the fit = ',1p2e16.7,a1) end  dy(i)',a1) 180 format(1x,i3,1p4e15.7,a1) 190 formatad LSQW.BAS 9/9/1986bn Weighted Least Squares Data Reduction using Gaussian Elimination4bx Tested with NBS Test PolynomialMb Robert H. Shermanb 1931 Camino Manzana, Los Alamos, New Mexico 87544b AH,OZ : INbMX2 : EF : MD : xbPR$" +#.#####^^^^" : PR1$"##" : PR2$" +#.#######^^^^"@cPR3$" +#.###^^^^" : PR4$" +#.######^^^^" : PR5$" +#.####^^^^"xc X(MX),Y(MX),XP(MX),WT(MX),DEV(MX),YC(MX),PCTD(MX)cWVOLD : SDOLDcQMD : A(Q,Q),R(Q),V(Q),AINV(Q,Q),SDV(Q)cQMD : P(Q)c ()"9/9/1986"/d "Weighted Least Squares Polynomial, Y= a + b*X + c*X^2 + ..."Md "(c) 1986 R.H. Sherman"d "Currently set to accept ";MX;" data pairs." : d "Create a Data File? Y/N =N ";QCR$d QCR$"Y" QCR$"y"  e" d : "*.DAT" : x : : "Input Date :";DT$e, T/e6J : "*.DAT" : fe@ "Input Data Filename FNM (.DAT) :",XFNM$ : eJ "i",#,XFNM$".DAT"eT #,HDR$ : HDR$e^ "Input Data Order : X,Y ( Y/N/Other=Y ) ";QD$fh "Are Weights included - (Y/N/Other=N) ";QW$ : QW$"y" QW$"Y"-fr QD$"n" QD$"N"ff| "Do you wish to scale the X's ( Y/N =N )";QSX$f "Do you wish to scale the Y's ( Y/N =N )";QSY$f QSX$"y" QSX$"Y"f QSY$"y" QSY$"Y"g QSX$"Y" "X Scale Factor =",SCFTRX7g QSY$"Y" "Y Scale Factor (Y=Y/scf)=",SCFTRYAgJJVg QW$ "Y" g QD$"N" #,X(J),Y(J) : #,Y(J),X(J)gWT(J) : g QD$"N" #,X(J),Y(J),WT(J) : #,Y(J),X(J),WT(J)g X(J)EF JJ : &0h JMX "** No More Data Allowed **" : &Rh QSX$"Y" X(J)X(J)SCFTRXvh QSY$"Y" Y(J)Y(J)SCFTRYh () &h h&NPJ h0 #h: NP "**FATAL ERROR** No Data Entered" : hD NP;" data points entered")iN "Degree of polynomial to be fitted ";D _iX D "**ERROR** Degree must be => 0" : NzibD(D) : DNP vil "**ERROR** Not enough data": NivD2Di DMD "**ERROR** degree too high": NiNDj I NP : XP(I) : I1j I D2 : I N k*  );HDR$;" ";DT$ : k4 :  );"X POWER COEFFICIENT SD OF COEFFICIENT" : "l> J N : " ";J, : PR2$;V(J); : ); : PR2$;SDV(J) : : 9lH QP$"Y" fGlR J Nl\  );J;) : PR2$;V(J); : ); : PR2$;SDV(J) : : lf J NPlp PR1$;J; : PR$;(WT(J));(X(J)),(Y(J)),(YC(J)),lz PR3$;(DEV(J)) m QP$"Y" Om  ); : PR1$;J; : PR4$;WT(J);X(J),Y(J),YC(J),DEV(J)Umm : "Weighted Variance, SD ="; : PR5$;WVAR,SDm ); : PR5$;WVOLD,SDOLDmSDOLDSD : WVOLDWVARm QP$"Y" n :  );"Weighted Variance, SD ="; : PR5$;WVAR,SD+n : : Pn : "Hit To Continue";QC$qn : "Continuation Options"n" 1 - Determine Specific Points"n" 2 - Fit another polynomial to the same data"n" 3 - Fit another set of data"o" 4 - TERMINATE program",o" 5 - Unity Weights"Oo$" 6 - Weights as 1/(Dev)**2"qo." 7 - Set Specific Weights"o8" 8 - Plot Errors"oB" 81 - Plot Percent Errors"oL" 9 - Reprint Data"oV" 10 - Output Fitted Parameters"p`" 11 - Restart Program"$pj" 99 - Return to CP/M"Qpt"What Next";Q : Q(Q) : Q cp~ Q rup Q p Q p Q QQ p Q  p Q p Q p Q ,q Q T : Nq Q !q Qc 2q Q  for YES ";QQ$q  QQ$""  : hq : "Enter ";EF;" to leave this mode"q : "X=";XV : XVEF q(YV : KN  ,r2YVYVXVV(K) : : "y= "; : PR2$;YV6r< erF Matrix Inversion by Gaussian EliminationorPLDrZ I L : J L : IJ B(I,J) : B(I,J)rd J,Irn : "Coefficient Matrix" sx I L : J L : (A(I,J)), : J : : I.s : I L : I,R(I) : I\s J L : IJ L : A(I,J) ds I|s "Singular Matrix"s Js K L : SA(J,K) : A(J,K)A(I,K) : A(I,K)SsSB(J,K) : B(J,K)B(I,K) : B(I,K)S : KsTA(J,J)/t K L : A(J,K)TA(J,K) : B(J,K)TB(J,K) : KMt L1 L : L1J htTA(L1,J) : K LtA(L1,K)A(L1,K)TA(J,K)tB(L1,K)B(L1,K)TB(J,K) : Kt L1,Jt : "|A|'"t I L : J L : (B(I,J)), : J : : Iu"Ku, I L : V(I) : J L : V(I)V(I)B(I,J)R(J) : J : Iju6 I L : I,V(I) : Iu@ "Press any character to continue ...",CH$uJuT "Printed Output ( Y/N =N )";QP$ : QP$"y" QP$"Y"u^uh vr Set all weights to unity/v| I NP : WT(I) : I : Rv Set weights to inverse of SDzv I NP : WT(I)DEV(I) : IvWMAXWT() : I NP : WT(I)WMAX WMAXWT(I)v Iv WMAXH w I NP : WT(I)WT(I)HWMAX : I : w Input individual weightsNw "Input Index of Weight to be set : ";IWw "Input Value for weight, W(";IW;") = "; : WT(IW) : w Create Least Squares Data File for Program LSQw R.H. Sherman 5/14/84)x : "Create a WEIGHTED Data File for Least Squares Regression Analysis"ex  ); : "Input File Name Xfnm(.DAT) ";XFNM$ : x "O",#,XFNM$".DAT"x "HEADER : ";HDR$ : #,HDR$x "Input Data as X,Y,Wt triplets"x& "Input 999,999 when done"x0I"y: "i= ";I; : " x,y,w = ",X(I),Y(I),WT(I)DyD X(I) Y(I) XVyNII : :yX N I : N;);X(N););Y(N););WT(N) : Nyb "Data OK? (Y/N) =Y ";Q$ : Q$"Y" Q$"y" Q$"" zl "Index of Data to be Corrected, i=",J:zv" Input Data as X,Y,Wt " : X(J),Y(J),WT(J) : Xhz N I : #, X(N),Y(N),WT(N) : Nqz #z "Finished with Data Input and File ( ";XFNM$;".DAT ) Written to Disk"zzzLONG(XFNM$) : T$""{ I LONG : C$(XFNM$,I,)3{ C$"a" C$"z" C$((C$) )U{T$T$C$ : I : XFNM$T$ : v{ ())"ERROR PLOT" : {EMINOL : EMAX EMIN :INITIAL VALUES{ I NP : Q PCTD(I)DEV(I) : PCTD(I)dDEV(I)Y(I) | PCTD(I)EMAX EMAXPCTD(I),| PCTD(I)EMIN EMINPCTD(I)4|  I| (EMIN)EMAX EMAX(EMIN) : EMINEMAX : MAKE ERROR LIMITS SYMMETRICAL FOR GRAPHING| MASK$"* | *"|*F$"###. +#.#######^^^^ +#.#######^^^^ +#.#####^^^^"X}4M1$((EMIN)):M2$((EMAX)):L(M1$)(M2$):LL: L(MASK$) L(MASK$)L : L}>M1$( (M1$))M1$(L)M2$: M1$ :MIN AND MAX VALUES}H  ) ((MASK$),"*")}R I NP~\M$MASK$:(M$,((PCTD(I)EMIN)((MASK$))(EMAXEMIN)),)"+"~f  ) M$%~p IA~z  ) ((MASK$),"*")O~ M1$ : v~ "Hit any key to continue ... ";~PL$ : PL$""  : ~ Output Fitted Parameters~ HDR$~ XFN$ "Input Filename for PARAMETER OUTPUT FNM (.PAR) :",XFNM$ : , "O",#,XFNM$".PAR": #,HDR$w I N : #, PR2$;V(I) : I,V(I) : I : # :  Compute Standard Deviation of fit and of fitted parameters J NP : YC(J) : JQ : J NP : QQY(J) : : XMQNP : T : GV J NP : YC(J) : KN  : YC(J)YC(J)X(J)V(K) : l DEV(J)Y(J)YC(J) TTWT(J)(DEV(J)) : GG(Y(J)XM) : $ WVART(NPD)Ȁ. G Td : L ؀8 SD(WVAR)B I N : SDV(I)SD(B(I,I)) : I L Y(J)XM) : $ WVART(NPD)Date: 9/11/1986 This library contains the source code for two least squares programs for fitting polynomials to data. There are some serious defects in most of the similar programs previously available in the library. The solution of the least squares problem requires the inversion of a matrix which is generally ill-conditioned. There are many terms whose absolute values are large, but whose sums are nearly zero. Therefore it is vitally important that every precaution be taken to see that the simple arithmetic is performed in a manner to preserve as much accuracy as is possible. For example, X(i)^n is usually performed by exp(n*log(X(i))) which is not as accurate as X(i)*X(i)*... since we are dealing with integer values of n. The codes presented here use double precision arithmetic. The programs may be tested by means of data generated from known polynomials of the type suggested several years ago in a paper from the National Bureau of Standards. Two such sets of data are included. NBS1.DAT is a set of numbers generated from a simple polynomial together with weights (unity). For data set RNDNRM5.DAT each value of Y(i) was perturbed by the addition of a value taken from a set of random normal deviates with mean of zero and a standard distribution of 1 ( and divided by 10000 to give s.d.=1.e-4). If a least squares code performs well it should reproduce the coefficients of the polynomial from which the data were generated with high accuracy. In the case of RNDNRM5.DAT the "standard deviation of the fit" should be close to the value of the built in "error". Both codes contained here pass these tests. LSQW.BAS is designed on a core culled from various other codes, but with modifications for accuracy and incorporating features which I find useful in my work. The console messages should be self explanitory. Code is included to permit the construction of a data file. If generated by an editor the format is: 1) a single title line 2) multiple data lines in the form X(i),Y(i),Weight(i) The code will sense an end-of-file condition and count the number of data items in the file. LSQF.FOR and LSQF.COM are the source and object codes for a somewhat faster and slightly more accurate least squares polynomial fit, but without many of the bells and whistles of LSQW.BAS. The format for the data file is: 2F12.0, preceeded by a single title line (10A8). The code will detect and end-of-filecondition for counting the number of data entries. If one chooses to use to use keyboard entry, the format is F12.0 ( which will permit properly justified exponential format. A file name may be in either case as conversion to upper case is included. Filename entry is a string of 11 legal filename characters as follows: 12345678901 filenameext ---> filename.ext It is hoped that you will find these programs useful. If you have any comments, corrections, upgrades, etc. I would be pleased to hear from you. Robert H. Sherman 1931 Camino Manzana Los Alamos, New Mexico 87544 . I would be pleased to hear from you. Robert H. Sherman 1931 NBS Polynomial w/ Random Normal Deviations, sd=1.e-4, # coef. = 6 (n=5) 1. 1001.110967 1.0 2. 1002.499328 1.0 3. 1004.275362 1.0 4. 1006.598470 1.0 5. 1009.687510 1.0 6. 1013.833462 1.0 7. 1019.411643 1.0 8. 1026.892810 1.0 9. 1036.856039 1.0 10. 1050.000125 1.0 11. 1067.156176 1.0 12. 1089.299114 1.0 13. 1117.560222 1.0 14. 1153.238404 1.0 15. 1197.812762 1.0 16. 1252.953558 1.0 17. 1320.536805 1.0 18. 1402.652900 1.0 19. 1501.620917 1.0 20. 1620.000039 1.0 21. 1760.601070 1.0 22. 1926.499168 1.0 23. 2121.045490 1.0 24. 2347.878322 1.0 25. 2610.937562 1.0 26. 2914.473447 1.0 27. 3263.061727 1.0 28. 3661.612898 1.0 29. 4115.385777 1.0 30. 4629.999974 1.0 31. 5211.446124 1.0 32. 5866.099171 1.0 33. 6600.730250 1.0 34. 7422.518440 1.0 35. 8339.062440 1.0 36. 9358.393736 1.0 37. 10488.986748 1.0 38. 11739.772797 1.0 39. 13120.150752 1.0 40. 14639.999874 1.0  10488.986748 1.0 38. 11739.772797 1.0 39. 13120.150752 1.0 40. NBS Polynomial with 4 Coefficients, n=3 (weights no included) 1. 0.0 2. -5.0 3. -20. 4. -51. 5. -104. 6. -185. 7. -300. 8. -455. 9. -656. 10. -909.