my.nls <- function(formula, data = sys.parent(), 
		   start = if(!is.null(attr(data, "parameters"))) 
		   attr(data, "parameters") else nls.initial(), control, 
		   algorithm = "default", trace = F,returnZ = F)
{
  convert.twiddle <- function(formula)
    {
      if(length(formula) < 3)
	return(formula)
      form <- call("~", call("-", formula[[2]], formula[[3]]))
      attr(form, "class") <- "formula"
      form
    }
  if(is.numeric(data))
    data <- sys.frame(data)
  cl <- class(data)
  if(inherits(data, "pframe")) {
    class(data) <- NULL
    pp <- parameters(data)
    if(length(pp)) {
      np <- names(pp)
      if(any(match(np, names(data), 0)))
	stop(
	     "can't have variables, parameters with same name"
	     )
			data[np] <- pp
    }
  }
    else if(inherits(data, "data.frame"))
      cl <- c("pframe", cl)
  class(data) <- NULL	
  ## First, figure out the par. names, make start a list
  switch(mode(start),
	 numeric = {
	   .parameters <- names(start)
	   start <- as.list(start)
	   names(start) <- .parameters
	 }
	 ,
	 list = {
	   .parameters <- names(start)
	 }
	 ,
	 NULL = .parameters <- parameter.names(formula, data),
	 stop("\"start\" should be numeric or list"))
  if(!length(.parameters))
    stop("names for parameters needed, from formula or from start")
  pn <- .parameters
  .expr <- formula	
  ## select the algorithm and possibly massage the formula
  if(length(start)) data[.parameters] <- start	#used by setup_nonlin
  data$.parameters <- .parameters
  nl.frame <- new.frame(data, F)
  frame.attr("class", nl.frame) <- cl	# in case data is returned
  formula <- switch(algorithm,
		    plinear = {
			if(length(formula) < 3)
			  stop("formula for plinear algorithm be of the form resp ~ array" )
			response <- eval(formula[[2]], nl.frame)
			design <- eval(formula[[3]], nl.frame)
			nnobs <- length(response)
			nlinear <- if(is.matrix(design)) dim(design)[2] else 
			length(design)/nnobs
			form <- call("~", formula[[3]])
			attr(form, "class") <- "formula"
			form
		      }
		    ,
		    default = convert.twiddle(formula))
  dims <- .C("setup_nonlin",
	     n = integer(3),
	     list(formula),
	     as.integer(nl.frame))$n
  npar <- dims[1]
  nderiv <- dims[2]
  nobs <- dims[3]
  resp <- eval(.expr[[2]], nl.frame)
  if(length(.expr) == 2) resp[] <- 0	
  ## setup_nonlin will have set .parameters if missing
  if(is.null(start)) {
    start <- list()
    if(is.null(names(start)) && length(start) == length(pn))
      names(start) <- pn
    for(i in .parameters)
      start[[i]] <- get(i, frame = nl.frame)
  }
  asgn <- start
  si <- 1
  for(i in 1:length(asgn)) {
    ni <- length(asgn[[i]])
    asgn[[i]] <- seq(from = si, length = ni)
    si <- si + ni
  }
  start <- unlist(start)
  controlvals <- nls.control()
  if(!missing(control))
    controlvals[names(control)] <- control
  ret.data <- is.logical(ret.data <- controlvals$data) && ret.data
  max.iterations <- if(is.null(controlvals$maxiter)) 50 * npar else 
  controlvals$maxiter
  settings <- c(0, max.iterations, controlvals$minscale, controlvals$
		tolerance, 0)
  if(algorithm == "plinear") {
    settings[1] <- 1
		dims <- c(dims, nlinear)
    start <- c(start, numeric(nlinear))
    pn <- c(pn, paste(".lin", 1:nlinear, sep = ""))
    dims[3] <- nnobs
		outmat <- array(c(response, numeric(nnobs * (npar + nlinear))), 
				c(nnobs, npar + nlinear + 1))
  }
    else outmat <- array(0, c(nobs, npar + 1))
  storage.mode(outmat) <- "double"
  storage.mode(start) <- "double"
  nls.trace <- if(missing(trace)) controlvals$trace else trace
  std.trace <- FALSE
  if(is.logical(nls.trace)) {
    if(std.trace <- nls.trace)
      nls.trace <- "trace.nls"
      else nls.trace <- NULL
  }
    else std.trace <- is.character(nls.trace) && nls.trace == "trace.nls"
  if(std.trace) {
    assign("trace.mat", array(0, c(max.iterations, npar + 2), list(
								   NULL, c("obj.", "conv.", paste("par", 1:npar)))), frame
	   = 1)
    assign("trace.expr", expression(trace.mat[last.iteration,  ] <- 
	it.row), frame = 1)
	}
  z <- .C("do_nls",
	  parameters = start,
	  dims = as.integer(dims),
	  control = as.double(settings),
	  outmat = outmat,
	  trace = list(nls.trace))
  if(z$control[5]){
    cat(switch(as.integer(z$control[5]),
	       "step factor reduced below minimum",
	       "maximum number of iterations exceeded",
	       "singular gradient matrix",
	       "singular design matrix",
	       "singular gradient matrix"))
    BADFLAG <- T}
  else
    BADFLAG <- F
  ##I added the BADFLAG thing so that instead of giving an error it would just
  ##send a warning. This is useful when you are calling the function from a
  ##loop. Also I added the gradient matrix to the things returned.
nls.out <- list(parameters = z$parameters, formula = .expr, call = 
		match.call(), residuals = z$outmat[, 1],
		Z=z$outmat[, -1],badflag=BADFLAG)
  if(algorithm == "plinear") {
    class(nls.out) <- c("nls.pl", "nls")
    R <- qr(z$outmat[, -1])$qr[1:(npar + nlinear),  , drop = F]
  }
    else {
      class(nls.out) <- "nls"
      R <- qr(z$outmat[, -1])$qr[1:npar,  , drop = F]
    }
  R[lower.tri(R)] <- 0
  nls.out$R <- R
  nls.out$fitted.values <- resp - nls.out$residuals
  nls.out$assign <- asgn
  if(ret.data) {
    data <- sys.frame(nl.frame)
    data$.parameters <- NULL	
###dbdetach does the right thing--only matters that nl.frame>1
    if(inherits(data, "pframe"))
      data <- dbdetach(data, nl.frame)
    nls.out$data <- data
  }
  if(std.trace && exists("last.iteration", frame = 1))
    nls.out$trace <- get("trace.mat", frame = 1)[1:get(
					"last.iteration", frame = 1),  ]
  nls.out
}