00001 
00002 
00003 
00004 
00005 
00006 #include "f2c.h"
00007 
00008 
00009 
00010 struct {
00011     doublereal xx[32768], yy[32768];
00012     integer numpts;
00013 } xydata_;
00014 
00015 #define xydata_1 xydata_
00016 
00017 struct {
00018     char ccode[4096];
00019     integer numcod, ivevar, numpar, ifxvar[26];
00020 } exdata_;
00021 
00022 #define exdata_1 exdata_
00023 
00024 
00025 
00026 static integer c__1 = 1;
00027 static integer c__5 = 5;
00028 static logical c_true = TRUE_;
00029 static integer c__2 = 2;
00030 static integer c_b73 = 983040;
00031 static integer c__32768 = 32768;
00032 
00033  MAIN__(void)
00034 {
00035     
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     
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     
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      int s_copy(char *, char *, ftnlen, ftnlen), s_cat(char *,
00080              char **, integer *, integer *, ftnlen);
00081     double sqrt(doublereal);
00082 
00083     
00084     static integer iend;
00085     static doublereal fvec[32768];
00086     static integer info;
00087     extern  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]      ;
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  int dnls1e_(U_fp, integer *, integer *, integer *,
00100              doublereal *, doublereal *, doublereal *, integer *, integer *, 
00101             integer *, doublereal *, integer *);
00102     extern  int ff_();
00103     static logical lchuse[26];
00104     extern  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     
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 
00140 
00141 
00142 
00143 
00144 
00145 
00146 
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 
00156 
00157 
00158 
00159 
00160 
00161 
00162 
00163 
00164 
00165 
00166 
00167 
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 
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 
00231     }
00232 
00233     if (iend <= istart) {
00234         iend = 32768;
00235     } else {
00236 
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 
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 
00309 
00310     for (i__ = 1; i__ <= 26; ++i__) {
00311         lchuse[i__ - 1] = FALSE_;
00312 
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 
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 
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 
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 
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 
00376 
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 
00408     }
00409 
00410 
00411 
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 
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 
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 
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 
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 
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 } 
00522 
00523 
00524 
00525 
00526  int ff_(integer *iflag, integer *m, integer *n, doublereal *
00527         x, doublereal *fvec, doublereal *fjac, integer *ldfjac)
00528 {
00529     
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     
00534     integer x_dim1, i__1, i__2;
00535 
00536     
00537     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00538      int s_stop(char *, ftnlen);
00539 
00540     
00541     static integer i__, j, k;
00542     static doublereal vaz[851968]       ;
00543     extern  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     
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 
00559 
00560 
00561 
00562 
00563 
00564 
00565 
00566 
00567 
00568 
00569 
00570 
00571 
00572 
00573 
00574 
00575 
00576     
00577     --fvec;
00578     x_dim1 = *n;
00579     --x;
00580 
00581     
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 
00596             }
00597 
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 
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 
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 } 
00627 
00628  int colfit_ () { MAIN__ (); return 0; }