aboutsummaryrefslogtreecommitdiff
path: root/math/nlfit/nliter.gx
diff options
context:
space:
mode:
authorJoe Hunkeler <jhunkeler@gmail.com>2015-08-11 16:51:37 -0400
committerJoe Hunkeler <jhunkeler@gmail.com>2015-08-11 16:51:37 -0400
commit40e5a5811c6ffce9b0974e93cdd927cbcf60c157 (patch)
tree4464880c571602d54f6ae114729bf62a89518057 /math/nlfit/nliter.gx
downloadiraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'math/nlfit/nliter.gx')
-rw-r--r--math/nlfit/nliter.gx59
1 files changed, 59 insertions, 0 deletions
diff --git a/math/nlfit/nliter.gx b/math/nlfit/nliter.gx
new file mode 100644
index 00000000..7857f852
--- /dev/null
+++ b/math/nlfit/nliter.gx
@@ -0,0 +1,59 @@
+include <mach.h>
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLITER - Routine to compute one iteration of the fitting process
+
+procedure nliter$t (nl, x, z, w, npts, nvars, ier)
+
+pointer nl # pointer to nl fitting structure
+PIXEL x[ARB] # independent variables (npts * nvars)
+PIXEL z[ARB] # function values (npts)
+PIXEL w[ARB] # weights (npts)
+int npts # number of points
+int nvars # number of independent variables
+int ier # error code
+
+int i, index
+PIXEL nlacpts$t(), nlresid$t()
+
+begin
+ # Do the initial accumulation.
+ NL_OLDSQ(nl) = nlacpts$t (nl, x, z, w, npts, nvars)
+
+ # Set up temporary parameter array.
+ call amov$t (PARAM(NL_PARAM(nl)), TRY(NL_TRY(nl)), NL_NPARAMS(nl))
+
+ repeat {
+
+ # Solve the matrix.
+ call nlsolve$t (nl, ier)
+ if (ier == NO_DEG_FREEDOM)
+ return
+
+ # Increment the parameters.
+ do i = 1, NL_NFPARAMS(nl) {
+ index = PLIST(NL_PLIST(nl)+i-1)
+ TRY(NL_TRY(nl)+index-1) = TRY(NL_TRY(nl)+index-1) +
+ DPARAM(NL_DPARAM(nl)+i-1)
+ }
+
+ # Reaccumulate the residuals and increment lambda.
+ NL_SUMSQ(nl) = nlresid$t (nl, x, z, w, npts, nvars)
+ #if (NL_OLDSQ(nl) > (NL_SUMSQ(nl) + NL_TOL(nl) * NL_SUMSQ(nl))) {
+ if (NL_OLDSQ(nl) >= NL_SUMSQ(nl)) {
+ call amov$t (TRY(NL_TRY(nl)), PARAM(NL_PARAM(nl)),
+ NL_NPARAMS(nl))
+ NL_LAMBDA(nl) = PIXEL (0.10) * NL_LAMBDA(nl)
+ break
+ } else
+ NL_LAMBDA(nl) = min (PIXEL(LAMBDAMAX),
+ PIXEL (10.0) * NL_LAMBDA(nl))
+
+ } until (NL_LAMBDA(nl) <= PIXEL(0.0) ||
+ NL_LAMBDA(nl) >= PIXEL(LAMBDAMAX))
+end