00628 { MAIN__ (); return 0; }

00528 {
00529     /* Format strings */
00530     static char fmt_101[] = "(\002 --- FF X=\002,5(1x,1pg12.5))";
00531     static char fmt_999[] = "(\002   *** FF has illegal IFLAG=\002,i5)";
00532 
00533     /* System generated locals */
00534     integer x_dim1, i__1, i__2;
00535 
00536     /* Builtin functions */
00537     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00538     /* Subroutine */ int s_stop(char *, ftnlen);
00539 
00540     /* Local variables */
00541     static integer i__, j, k;
00542     static doublereal vaz[851968]       /* was [32768][26] */;
00543     extern /* Subroutine */ int parevec_(integer *, char *, doublereal *, 
00544             doublereal *, doublereal *, doublereal *, doublereal *, 
00545             doublereal *, doublereal *, doublereal *, doublereal *, 
00546             doublereal *, doublereal *, doublereal *, doublereal *, 
00547             doublereal *, doublereal *, doublereal *, doublereal *, 
00548             doublereal *, doublereal *, doublereal *, doublereal *, 
00549             doublereal *, doublereal *, doublereal *, doublereal *, 
00550             doublereal *, integer *, doublereal *, ftnlen);
00551 
00552     /* Fortran I/O blocks */
00553     static cilist io___54 = { 0, 6, 0, fmt_101, 0 };
00554     static cilist io___59 = { 0, 6, 0, fmt_999, 0 };
00555 
00556 
00557 
00558 /*  Routine supplied to DNLS1E to compute the functions we are */
00559 /*  are trying to fit. */
00560 
00561 /* .......................................................................
00562  */
00563 /*  Common block XYDATA holds the data read from disk */
00564 
00565 /* .......................................................................
00566  */
00567 /*  Common block EXDATA holds the parsed expression for evaluation */
00568 
00569 /* .......................................................................
00570  */
00571 /*  Local variables */
00572 
00573 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00574  */
00575 
00576     /* Parameter adjustments */
00577     --fvec;
00578     x_dim1 = *n;
00579     --x;
00580 
00581     /* Function Body */
00582     if (*iflag == 0) {
00583         s_wsfe(&io___54);
00584         do_fio(&x_dim1, (char *)&x[1], (ftnlen)sizeof(doublereal));
00585         e_wsfe();
00586 /* ...................................................................
00587 .... */
00588     } else if (*iflag == 1) {
00589         i__1 = *n;
00590         for (i__ = 1; i__ <= i__1; ++i__) {
00591             k = exdata_1.ifxvar[i__ - 1];
00592             i__2 = *m;
00593             for (j = 1; j <= i__2; ++j) {
00594                 vaz[j + (k << 15) - 32769] = x[i__];
00595 /* L105: */
00596             }
00597 /* L110: */
00598         }
00599         i__1 = *m;
00600         for (j = 1; j <= i__1; ++j) {
00601             vaz[j + (exdata_1.ivevar << 15) - 32769] = xydata_1.xx[j - 1];
00602 /* L115: */
00603         }
00604         parevec_(&exdata_1.numcod, exdata_1.ccode, vaz, &vaz[32768], &vaz[
00605                 65536], &vaz[98304], &vaz[131072], &vaz[163840], &vaz[196608],
00606                  &vaz[229376], &vaz[262144], &vaz[294912], &vaz[327680], &vaz[
00607                 360448], &vaz[393216], &vaz[425984], &vaz[458752], &vaz[
00608                 491520], &vaz[524288], &vaz[557056], &vaz[589824], &vaz[
00609                 622592], &vaz[655360], &vaz[688128], &vaz[720896], &vaz[
00610                 753664], &vaz[786432], &vaz[819200], m, &fvec[1], 8L);
00611         i__1 = *m;
00612         for (i__ = 1; i__ <= i__1; ++i__) {
00613             fvec[i__] -= xydata_1.yy[i__ - 1];
00614 /* L120: */
00615         }
00616 /* ...................................................................
00617 .... */
00618     } else {
00619         s_wsfe(&io___59);
00620         do_fio(&c__1, (char *)&(*iflag), (ftnlen)sizeof(integer));
00621         e_wsfe();
00622         s_stop("", 0L);
00623     }
00624 
00625     return 0;
00626 } /* ff_ */

00034 {
00035     /* Format strings */
00036     static char fmt_101[] = "(\002 Enter data file name : \002$)";
00037     static char fmt_111[] = "(a)";
00038     static char fmt_102[] = "(\002 Start, stop indices in data file  [0,0=al"
00039             "l] ? \002$)";
00040     static char fmt_112[] = "(2i10)";
00041     static char fmt_119[] = "(\002   *** failure when skipping record \002,i"
00042             "5)";
00043     static char fmt_201[] = "(\002   --- read\002,i6,\002 data lines from fi"
00044             "le\002)";
00045     static char fmt_301[] = "(/\002 Expression to fit 2nd column to 1st ?"
00046             " \002$)";
00047     static char fmt_309[] = "(\002   *** expression too complex!\002)";
00048     static char fmt_338[] = "(\002   *** not enough variable names in expres"
00049             "sion!\002)";
00050     static char fmt_339[] = "(\002   *** too many variable names in expressi"
00051             "on!\002)";
00052     static char fmt_341[] = "(\002 Which variable represents 1st data column"
00053             " ? \002/1x,a,\002 : \002$)";
00054     static char fmt_351[] = "(\002  Initial value for parameter \002,a,\002 "
00055             "? \002$)";
00056     static char fmt_501[] = "(\002   --- \002,a,\002 exits with INFO =\002,i"
00057             "2)";
00058     static char fmt_511[] = "(3x,a,\002 -> \002,1pg14.7,\002 +/- \002,1pg14."
00059             "7)";
00060     static char fmt_521[] = "(\002   --- mean absolute deviation =\002,1pg10"
00061             ".3)";
00062     static char fmt_601[] = "(\002 Filename to save fit error curve in (blan"
00063             "k=none) ? \002$)";
00064     static char fmt_611[] = "(1pg20.13,1x,1pg20.13)";
00065 
00066     /* System generated locals */
00067     address a__1[2];
00068     integer i__1, i__2, i__3[2];
00069     doublereal d__1;
00070     char ch__1[1];
00071     olist o__1;
00072     cllist cl__1;
00073 
00074     /* Builtin functions */
00075     integer s_wsfe(cilist *), e_wsfe(void), s_rsfe(cilist *), do_fio(integer *
00076             , char *, ftnlen), e_rsfe(void), f_open(olist *), f_clos(cllist *)
00077             , s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
00078             e_rsle(void), s_cmp(char *, char *, ftnlen, ftnlen);
00079     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen), s_cat(char *,
00080              char **, integer *, integer *, ftnlen);
00081     double sqrt(doublereal);
00082 
00083     /* Local variables */
00084     static integer iend;
00085     static doublereal fvec[32768];
00086     static integer info;
00087     extern /* Subroutine */ int dcov_(U_fp, integer *, integer *, integer *, 
00088             doublereal *, doublereal *, doublereal *, integer *, integer *, 
00089             doublereal *, doublereal *, doublereal *, doublereal *);
00090     static integer ierr;
00091     static doublereal xpar[66];
00092     static integer iopt;
00093     static doublereal work[983040]      /* was [32768][30] */;
00094     static integer i__, m, n;
00095     static char cfile[64], chans[1];
00096     static integer numch;
00097     static char cexpr[128];
00098     static integer iwork[66];
00099     extern /* Subroutine */ int dnls1e_(U_fp, integer *, integer *, integer *,
00100              doublereal *, doublereal *, doublereal *, integer *, integer *, 
00101             integer *, doublereal *, integer *);
00102     extern /* Subroutine */ int ff_();
00103     static logical lchuse[26];
00104     extern /* Subroutine */ int parser_(char *, logical *, integer *, char *, 
00105             ftnlen, ftnlen);
00106     static integer ilower, iupper, istart, nprint, ich;
00107     static doublereal xin, yin, tol, war1[30], war2[30], war3[30], war4[32768]
00108             ;
00109 
00110     /* Fortran I/O blocks */
00111     static cilist io___1 = { 0, 6, 0, fmt_101, 0 };
00112     static cilist io___2 = { 0, 5, 0, fmt_111, 0 };
00113     static cilist io___5 = { 0, 6, 0, fmt_102, 0 };
00114     static cilist io___6 = { 0, 5, 0, fmt_112, 0 };
00115     static cilist io___10 = { 1, 1, 1, 0, 0 };
00116     static cilist io___13 = { 0, 6, 0, fmt_119, 0 };
00117     static cilist io___14 = { 1, 1, 1, 0, 0 };
00118     static cilist io___15 = { 0, 6, 0, fmt_201, 0 };
00119     static cilist io___16 = { 0, 6, 0, fmt_301, 0 };
00120     static cilist io___17 = { 1, 5, 1, fmt_111, 0 };
00121     static cilist io___19 = { 0, 6, 0, fmt_309, 0 };
00122     static cilist io___25 = { 0, 6, 0, fmt_338, 0 };
00123     static cilist io___26 = { 0, 6, 0, fmt_339, 0 };
00124     static cilist io___27 = { 0, 6, 0, fmt_341, 0 };
00125     static cilist io___28 = { 0, 5, 0, fmt_111, 0 };
00126     static cilist io___30 = { 0, 6, 0, fmt_351, 0 };
00127     static cilist io___31 = { 1, 5, 1, 0, 0 };
00128     static cilist io___42 = { 0, 6, 0, fmt_501, 0 };
00129     static cilist io___47 = { 0, 6, 0, fmt_501, 0 };
00130     static cilist io___48 = { 0, 6, 0, fmt_511, 0 };
00131     static cilist io___49 = { 0, 6, 0, fmt_511, 0 };
00132     static cilist io___50 = { 0, 6, 0, fmt_521, 0 };
00133     static cilist io___51 = { 0, 6, 0, fmt_601, 0 };
00134     static cilist io___52 = { 0, 5, 0, fmt_111, 0 };
00135     static cilist io___53 = { 0, 1, 0, fmt_611, 0 };
00136 
00137 
00138 
00139 /*  Program to do a nonlinear least squares fit of a two-column file in */
00140 /*  the format */
00141 
00142 /*    x1  y1 */
00143 /*    x2  y2 */
00144 /*    ..  .. */
00145 /*    xN  yN */
00146 
00147 /*  to a formula of the form  y(x) = expression(x,a,b,c,...), where */
00148 /*  a, b, c, ... are parameters to be found. */
00149 /* -----------------------------------------------------------------------
00150  */
00151 /*  Common block XYDATA holds the data read from disk */
00152 
00153 /* .......................................................................
00154  */
00155 /*  Common block EXDATA holds the parsed expression for evaluation */
00156 
00157 /* .......................................................................
00158  */
00159 /*  local variables */
00160 
00161 
00162 
00163 
00164 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00165  */
00166 
00167 /*  Open data file */
00168 
00169 L100:
00170     s_wsfe(&io___1);
00171     e_wsfe();
00172     s_rsfe(&io___2);
00173     do_fio(&c__1, cfile, 64L);
00174     e_rsfe();
00175     if (*(unsigned char *)cfile == ' ') {
00176         goto L9900;
00177     }
00178 
00179     o__1.oerr = 1;
00180     o__1.ounit = 1;
00181     o__1.ofnmlen = 64;
00182     o__1.ofnm = cfile;
00183     o__1.orl = 0;
00184     o__1.osta = "OLD";
00185     o__1.oacc = 0;
00186     o__1.ofm = "FORMATTED";
00187     o__1.oblnk = 0;
00188     ierr = f_open(&o__1);
00189 
00190     if (ierr != 0) {
00191         cl__1.cerr = 0;
00192         cl__1.cunit = 1;
00193         cl__1.csta = 0;
00194         f_clos(&cl__1);
00195         goto L100;
00196     }
00197 /* .......................................................................
00198  */
00199 /*  Read data in */
00200 
00201     s_wsfe(&io___5);
00202     e_wsfe();
00203     s_rsfe(&io___6);
00204     do_fio(&c__1, (char *)&istart, (ftnlen)sizeof(integer));
00205     do_fio(&c__1, (char *)&iend, (ftnlen)sizeof(integer));
00206     e_rsfe();
00207 
00208     i__1 = istart - 1;
00209     for (i__ = 1; i__ <= i__1; ++i__) {
00210         ierr = s_rsle(&io___10);
00211         if (ierr != 0) {
00212             goto L100001;
00213         }
00214         ierr = do_lio(&c__5, &c__1, (char *)&xin, (ftnlen)sizeof(doublereal));
00215         if (ierr != 0) {
00216             goto L100001;
00217         }
00218         ierr = do_lio(&c__5, &c__1, (char *)&yin, (ftnlen)sizeof(doublereal));
00219         if (ierr != 0) {
00220             goto L100001;
00221         }
00222         ierr = e_rsle();
00223 L100001:
00224         if (ierr != 0) {
00225             s_wsfe(&io___13);
00226             do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
00227             e_wsfe();
00228             goto L9900;
00229         }
00230 /* L110: */
00231     }
00232 
00233     if (iend <= istart) {
00234         iend = 32768;
00235     } else {
00236 /* Computing MIN */
00237         i__1 = 32768, i__2 = iend - istart + 1;
00238         iend = min(i__1,i__2);
00239     }
00240 
00241     i__ = 0;
00242 L200:
00243     ierr = s_rsle(&io___14);
00244     if (ierr != 0) {
00245         goto L100002;
00246     }
00247     ierr = do_lio(&c__5, &c__1, (char *)&xin, (ftnlen)sizeof(doublereal));
00248     if (ierr != 0) {
00249         goto L100002;
00250     }
00251     ierr = do_lio(&c__5, &c__1, (char *)&yin, (ftnlen)sizeof(doublereal));
00252     if (ierr != 0) {
00253         goto L100002;
00254     }
00255     ierr = e_rsle();
00256 L100002:
00257     if (ierr == 0) {
00258         ++i__;
00259         xydata_1.xx[i__ - 1] = xin;
00260         xydata_1.yy[i__ - 1] = yin;
00261         if (i__ < iend) {
00262             goto L200;
00263         }
00264     }
00265 
00266     cl__1.cerr = 0;
00267     cl__1.cunit = 1;
00268     cl__1.csta = 0;
00269     f_clos(&cl__1);
00270     xydata_1.numpts = i__;
00271     s_wsfe(&io___15);
00272     do_fio(&c__1, (char *)&xydata_1.numpts, (ftnlen)sizeof(integer));
00273     e_wsfe();
00274     if (xydata_1.numpts <= 1) {
00275         goto L9900;
00276     }
00277 /* .......................................................................
00278  */
00279 /*  Read in expression that will be used to fit */
00280 
00281 L300:
00282     s_wsfe(&io___16);
00283     e_wsfe();
00284     ierr = s_rsfe(&io___17);
00285     if (ierr != 0) {
00286         goto L100003;
00287     }
00288     ierr = do_fio(&c__1, cexpr, 128L);
00289     if (ierr != 0) {
00290         goto L100003;
00291     }
00292     ierr = e_rsfe();
00293 L100003:
00294     if (ierr != 0) {
00295         goto L9900;
00296     }
00297     parser_(cexpr, &c_true, &exdata_1.numcod, exdata_1.ccode, 128L, 8L);
00298     if (exdata_1.numcod <= 0) {
00299         goto L300;
00300     }
00301     if (exdata_1.numcod > 512) {
00302         s_wsfe(&io___19);
00303         e_wsfe();
00304         goto L9900;
00305     }
00306 /* .......................................................................
00307  */
00308 /*  Get the names of all the variables referred to in the expression */
00309 
00310     for (i__ = 1; i__ <= 26; ++i__) {
00311         lchuse[i__ - 1] = FALSE_;
00312 /* L310: */
00313     }
00314 
00315     iupper = 'A' - 1;
00316     ilower = 'a' - 1;
00317     i__1 = exdata_1.numcod;
00318     for (i__ = 1; i__ <= i__1; ++i__) {
00319         if (s_cmp(exdata_1.ccode + (i__ - 1 << 3), "PUSHSYM", 8L, 7L) == 0) {
00320             ich = *(unsigned char *)&exdata_1.ccode[i__ * 8] - iupper;
00321             lchuse[ich - 1] = TRUE_;
00322         }
00323 /* L320: */
00324     }
00325 
00326     s_copy(cfile, " ", 64L, 1L);
00327     numch = 0;
00328     for (i__ = 1; i__ <= 26; ++i__) {
00329         if (lchuse[i__ - 1]) {
00330             ++numch;
00331             i__1 = (numch << 1) - 1;
00332 /* Writing concatenation */
00333             *(unsigned char *)&ch__1[0] = i__ + iupper;
00334             i__3[0] = 1, a__1[0] = ch__1;
00335             i__3[1] = 1, a__1[1] = " ";
00336             s_cat(cfile + i__1, a__1, i__3, &c__2, (numch << 1) + 1 - i__1);
00337         }
00338 /* L330: */
00339     }
00340 
00341     if (numch <= 1) {
00342         s_wsfe(&io___25);
00343         e_wsfe();
00344         goto L300;
00345     }
00346     if (numch > xydata_1.numpts) {
00347         s_wsfe(&io___26);
00348         e_wsfe();
00349         goto L300;
00350     }
00351 /* .......................................................................
00352  */
00353 /*  Find out which variable is the 1st data column */
00354 
00355 L340:
00356     s_wsfe(&io___27);
00357     do_fio(&c__1, cfile, numch << 1);
00358     e_wsfe();
00359     s_rsfe(&io___28);
00360     do_fio(&c__1, chans, 1L);
00361     e_rsfe();
00362     if (*(unsigned char *)chans >= 'A' && *(unsigned char *)chans <= 'Z') {
00363         ich = *(unsigned char *)chans - iupper;
00364     } else if (*(unsigned char *)chans >= 'a' && *(unsigned char *)chans <= 
00365             'z') {
00366         ich = *(unsigned char *)chans - ilower;
00367     } else {
00368         goto L340;
00369     }
00370     if (! lchuse[ich - 1]) {
00371         goto L340;
00372     }
00373 /* .......................................................................
00374  */
00375 /*  Assign 1st column variable name to VEctor VARiable, and */
00376 /*  all others to be FiXed Variables;  get initial values for the latter 
00377 */
00378 
00379     exdata_1.ivevar = ich;
00380     exdata_1.numpar = numch - 1;
00381     ich = 0;
00382     for (i__ = 1; i__ <= 26; ++i__) {
00383         if (lchuse[i__ - 1] && i__ != exdata_1.ivevar) {
00384             ++ich;
00385             exdata_1.ifxvar[ich - 1] = i__;
00386 
00387 L350:
00388             s_wsfe(&io___30);
00389             *(unsigned char *)&ch__1[0] = (char) (i__ + iupper);
00390             do_fio(&c__1, ch__1, 1L);
00391             e_wsfe();
00392             ierr = s_rsle(&io___31);
00393             if (ierr != 0) {
00394                 goto L100004;
00395             }
00396             ierr = do_lio(&c__5, &c__1, (char *)&xpar[ich - 1], (ftnlen)
00397                     sizeof(doublereal));
00398             if (ierr != 0) {
00399                 goto L100004;
00400             }
00401             ierr = e_rsle();
00402 L100004:
00403             if (ierr != 0) {
00404                 goto L350;
00405             }
00406         }
00407 /* L360: */
00408     }
00409 /* .......................................................................
00410  */
00411 /*  Call the fitting routine */
00412 
00413     iopt = 1;
00414     m = xydata_1.numpts;
00415     n = exdata_1.numpar;
00416     tol = 1e-6;
00417     nprint = 0;
00418 
00419 /*     ************ */
00420     dnls1e_((U_fp)ff_, &iopt, &m, &n, xpar, fvec, &tol, &nprint, &info, iwork,
00421              work, &c_b73);
00422 /*     ************ */
00423 
00424 /*  Write results out: */
00425 
00426     s_wsfe(&io___42);
00427     do_fio(&c__1, "DNLS1E", 6L);
00428     do_fio(&c__1, (char *)&info, (ftnlen)sizeof(integer));
00429     e_wsfe();
00430     if (info == 0) {
00431         goto L9900;
00432     }
00433 
00434     dcov_((U_fp)ff_, &iopt, &m, &n, xpar, fvec, work, &c__32768, &info, war1, 
00435             war2, war3, war4);
00436 
00437     s_wsfe(&io___47);
00438     do_fio(&c__1, "DCOV", 4L);
00439     do_fio(&c__1, (char *)&info, (ftnlen)sizeof(integer));
00440     e_wsfe();
00441     if (info == 0) {
00442         goto L9900;
00443     }
00444 
00445     i__1 = exdata_1.numpar;
00446     for (i__ = 1; i__ <= i__1; ++i__) {
00447         if (info == 1) {
00448             xin = sqrt(work[i__ + (i__ << 15) - 32769]);
00449             s_wsfe(&io___48);
00450             *(unsigned char *)&ch__1[0] = (char) (exdata_1.ifxvar[i__ - 1] + 
00451                     iupper);
00452             do_fio(&c__1, ch__1, 1L);
00453             do_fio(&c__1, (char *)&xpar[i__ - 1], (ftnlen)sizeof(doublereal));
00454             do_fio(&c__1, (char *)&xin, (ftnlen)sizeof(doublereal));
00455             e_wsfe();
00456         } else {
00457             s_wsfe(&io___49);
00458             *(unsigned char *)&ch__1[0] = (char) (exdata_1.ifxvar[i__ - 1] + 
00459                     iupper);
00460             do_fio(&c__1, ch__1, 1L);
00461             do_fio(&c__1, (char *)&xpar[i__ - 1], (ftnlen)sizeof(doublereal));
00462             e_wsfe();
00463         }
00464 /* L510: */
00465     }
00466 
00467     xin = 0.;
00468     i__1 = xydata_1.numpts;
00469     for (i__ = 1; i__ <= i__1; ++i__) {
00470         xin += (d__1 = fvec[i__ - 1], abs(d__1));
00471 /* L520: */
00472     }
00473     xin /= xydata_1.numpts;
00474     s_wsfe(&io___50);
00475     do_fio(&c__1, (char *)&xin, (ftnlen)sizeof(doublereal));
00476     e_wsfe();
00477 /* .......................................................................
00478  */
00479 /*  Write fit error curve to a file (if desired) */
00480 
00481 L600:
00482     s_wsfe(&io___51);
00483     e_wsfe();
00484     s_rsfe(&io___52);
00485     do_fio(&c__1, cfile, 64L);
00486     e_rsfe();
00487     if (*(unsigned char *)cfile != ' ') {
00488         o__1.oerr = 1;
00489         o__1.ounit = 1;
00490         o__1.ofnmlen = 64;
00491         o__1.ofnm = cfile;
00492         o__1.orl = 0;
00493         o__1.osta = "UNKNOWN";
00494         o__1.oacc = 0;
00495         o__1.ofm = "FORMATTED";
00496         o__1.oblnk = 0;
00497         ierr = f_open(&o__1);
00498         if (ierr != 0) {
00499             goto L600;
00500         }
00501 
00502         i__1 = xydata_1.numpts;
00503         for (i__ = 1; i__ <= i__1; ++i__) {
00504             s_wsfe(&io___53);
00505             do_fio(&c__1, (char *)&xydata_1.xx[i__ - 1], (ftnlen)sizeof(
00506                     doublereal));
00507             do_fio(&c__1, (char *)&fvec[i__ - 1], (ftnlen)sizeof(doublereal));
00508             e_wsfe();
00509 /* L610: */
00510         }
00511 
00512         cl__1.cerr = 0;
00513         cl__1.cunit = 1;
00514         cl__1.csta = 0;
00515         f_clos(&cl__1);
00516     }
00517 /* -----------------------------------------------------------------------
00518  */
00519 L9900:
00520     return 0;
00521 } /* MAIN__ */