.packageName <- "ftnonpar"
"frun"<-function(y,...,alpha=0.5,r=0,mr=0)
{
#	
# INPUTS:		
#	y:	the data
#    ... :	an optional argument which specifies
#		the approximate positions of the local extreme values. 
#		These should be consisten with the run length otherwise 
#		the result will be nonsensical. Should you wish to use
#		this option then you should first run the macro
#		without it. The item list$xb of the output list gives 
#		the acceptable limits of the local extreme values. You
#		can then specify the positions within these limits.
#   alpha:	Qunatile determining the acceptable run length.
#       r:	Acceptable run length:	 Overrides alpha if not 0
#      mr:	mr=0 minimizes the run length consistent with the
#		number of local extreme values found for the specified 
#		run length. mr=1 disables the option.
#
#		
	args<-list(y,...)
	k<-length(args)
	y<-args[[1]]
	nxl<-1
	il<-0
	if (k>1)
		{
		il<-args[[2]]
		nxl<-length(il)
		}
		
#	sample size n
	n<- length(y)
#
# fortran subroutine
#
    	tmp <- .Fortran(
		"frun", 
		as.double(y), 
		double(n), 
		double(n), 
		double(2*n),
		double(2*n),
		double(n),
		integer(n),
		integer(n),
		as.integer(n),
		as.integer(r),
		as.double(alpha),
		as.integer(mr),
		as.integer(nxl),
		integer(1),
                PACKAGE="ftnonpar"
		)
#
# OUTPUTS		
#
#	number of extremes
	nx<-tmp[[14]]
#
#	run length: may differ from specified if mr=1	 
	r<-tmp[[10]]
	
#    	lower bounds 
	l1<- tmp[[2]]
#
#       upper bound
	u1<- tmp[[3]]
#
#	bounds for location of extremes: the position of the ith
#	extreme value lies between xb[2*i-1] and xb[2*i]		 
	xb<-tmp[[7]]
	xb<-xb[1:(2*nx)]
#
#       lower bound with specified extremes: the default choices for
#	the positions of the local extreme values are the mid-points
#	of the intervals specified by xb above.			
	l2<- tmp[[4]]
	l2<-l2[1:n]
#
#	upper bound with specified extremes
 	u2<- tmp[[5]]
	u2<-u2[1:n]
#
#	function between l2 and u2 satisfying run condition
	f<- tmp[[6]]
#
#       IN GENERAL THE MEAN OF THE BOUNDS l2 AND u2 (l2+u2)/2
#	GIVES A BETTER REGRESSION FUNCTION THAN f. HOWEVER THIS
#	FUNCTION IS INFINITE AT THE TWO ENDPOINTS AND AT LOCAL EXTREME
#	VALUES. IN THESE INTERVALS IT CAN BE REPLACED BY ANY VALUES
#	WHICH DO NOT ALTER THE NUMBER OF LOCAL EXTREME VALUES. THE
#	MEDIAN OF THE y-VALUES IN THESE INTERVALS IS A REASONABLE 
#	CHOICE.
#			 				
	list(l1=l1,u1=u1,l2=l2,u2=u2,f=f,xb=xb,nx=nx,r=r)
#
}


"quantpmreg" <-
function(y,beta=0.5,squeezing.factor=0.5,verbose=FALSE,localsqueezing=TRUE,DYADIC=TRUE,thr.const=2,extrema.nr = -1,bandwidth= -1)
{
n<-length(y)
sigma <- 1

if (bandwidth < 0)
        firstlambda <- 2^(floor(log(length(y), base = 2)) - 1)
    else firstlambda <- bandwidth

lambda <- rep(firstlambda,n-1)
currprecision <- firstlambda

while(1<2)
  {
  tmp <- genstring(y,lambda,beta,1)
  y.string <- tmp$y
  y.mr <- genmrcheck(y-y.string,sigma=sigma,beta=beta,method=1,DYADIC=DYADIC,thr.const=thr.const)
  if(verbose)
    {
    par(mfrow=c(1,2))
    plot(y,col="lightgrey")
    lines(y.string,col="red")
    plot(lambda,ylim=c(-max(lambda),max(lambda)),ty="l")
    lines(-lambda)
    lines(cumsum(sign(y-y.string)),col="red")
    print(c("lambda=",min(lambda)))
    print("Press Enter")
    readline()
    }
  if (bandwidth > 0)
     break
  if (extrema.nr > 0) {
    if (tmp$kext > extrema.nr)
       lambda <- lambda + currprecision
            if (currprecision < 0.5) {
                if (tmp$kext <= extrema.nr)
                  break
            }
            else {
                currprecision <- currprecision/2
                lambda <- lambda - currprecision
            }
        }
        else {

  if(sum(y.mr)<0.1)
    break
  if(!localsqueezing)
    lambda <- lambda*squeezing.factor
  else
    lambda[y.mr[-1] | y.mr[-n]] <- lambda[y.mr[-1] | y.mr[-n]] * squeezing.factor
}

  }
list(y=y.string,lambda=lambda,nmax=tmp$kext)
}
"genstring" <-
function(y,lambda,beta=0.5,method=1)
{
n <- length(y)
if(length(lambda)==1)
  lambda <- c(rep(lambda,n-1),0)
else
if(length(lambda)==n-1)
  lambda <- c(lambda,0)
if(method == 1)
  {
  tmp <- sort(y)
  eps <- min(c(0.00001*(max(y)-min(y)),min(tmp[-1] - tmp[-n])/4))
  if(eps < 1e-36)
    {
    if(mad(y[-1]-y[-n])>0)
      y <- y + rnorm(n,0,0.001*mad(y[-1]-y[-n]))
    else
      y <- y + rnorm(n,0,1e-12)
    tmp <- sort(y)
    eps <- min(tmp[-1] - tmp[-n])/4
    }
  }
else
  eps <- 0

print(c("eps is",eps))
tmp <- .C("genstring",y=as.double(y),as.integer(length(y)),as.double(lambda),
as.double(beta),as.integer(method),as.double(eps),kext=as.integer(0),PACKAGE="ftnonpar")

list(y=tmp$y,kext=tmp$kext)
}
"genmrcheck" <-
function(res,thresh=-1,sigma=1,DYADIC=FALSE,beta=0.5,method=1,thr.const=2)
{
if(method == 1)
  {
  res[res> 1e-04] <- beta
  res[res< -1e-04] <- -1+beta
  sigma <- sqrt(beta*(1-beta))

  currq <- pnorm(sqrt(thr.const*log(length(res))))
  tmp <- ceiling(qbinom(1-currq,1:length(res),beta) )
  upperthresh <- -tmp+beta*(1:length(res))
  tmp <- floor(qbinom(currq,1:length(res),beta)) 
  lowerthresh <- -tmp+beta*(1:length(res))
  }
else
  {
  upperthresh <- sqrt(thr.const*log(length(res)))*sigma*sqrt(1:length(res))
  lowerthresh <- -sqrt(thr.const*log(length(res)))*sigma*sqrt(1:length(res))
  }

  .C("genmrcheck",res=as.double(c(0,cumsum(res))),as.integer(length(res)),as.double(lowerthresh),as.double(upperthresh),as.integer(DYADIC),PACKAGE="ftnonpar")$res[1:length(res)]
}


"l1pmreg" <-
function(y,squeezing.factor=0.5,verbose=FALSE,localsqueezing=TRUE,DYADIC=TRUE,thr.const=2,extrema.nr = -1,bandwidth= -1)
{
quantpmreg(y,0.5,squeezing.factor,verbose,localsqueezing,DYADIC,thr.const,extrema.nr,bandwidth)
}
"mintvmon" <-
function(y,sigma=-1,DYADIC=TRUE,thresh=-1,method=2,MONCONST=TRUE,CONVCONST=FALSE)
{
n <- length(y)
if(thresh<0)
  {
  if(sigma < 0)
    sigma <- mad(y[-1]-y[-n])/sqrt(2)
  print(sigma)
  thresh <- sqrt(2*log(n))*sigma
  }
tmp <- .C("mintvmon", as.integer(n),
f=double(n),derivsign=integer(n-1),secsign=integer(n-1),as.integer(method),as.integer(DYADIC),as.double(thresh),as.double(y),as.integer(MONCONST),as.integer(CONVCONST),jact=integer(n),kact=integer(n),signact=integer(n),nact=integer(1),piecesleft=integer(n),piecesright=integer(n),PACKAGE="ftnonpar")

list(y=tmp$f,derivsign=tmp$derivsign,secsign=tmp$secsign,jact=tmp$jact[1:tmp$nact],kact=tmp$kact[1:tmp$nact],signact=tmp$signact[1:tmp$nact],pl=tmp$piecesleft[1:tmp$nact],pr=tmp$piecesright[1:tmp$nact])

}
"pmden" <-
function (x, DISCR=FALSE,verbose = FALSE, bandwidth = -1, extrema.nr = -1, accuracy = mad(x)/1000, 
    extrema.mean = TRUE, maxkuipnr = 19, asympbounds = FALSE, tolerance = 0.001) 
{
    nsamp <- length(x)
    if (asympbounds || nsamp > max(kuipdiffbounds.x)) 
        currbounds <- grenzen[length(kuipdiffbounds.x), ] * sqrt(max(kuipdiffbounds.x))/sqrt(nsamp)
    else {
        currbounds <- double(maxkuipnr)
        for (i in 1:maxkuipnr) currbounds[i] <- approx(kuipdiffbounds.x, 
            kuipdiffbounds[, i], nsamp)$y
    }
    if (maxkuipnr > dim(kuipdiffbounds)[2]) 
        stop("maxkuipnr is too large")
    if(DISCR)
      {
      datax <- as.double(levels(as.factor(x)))
      N <- length(datax)
      dataemp <- c(0,as.double(summary(as.factor(x),maxsum=N))/nsamp)
      fdist.y <- cumsum(dataemp)
      x <- 0:N
      nsamp <- N+1
      }
    else
      {
      dataemp <- c(0,rep(1/(nsamp-1),nsamp-1))
      x <- sort(x)
      if (min((x[-1] - x[-nsamp])/(x[nsamp] - x[1])) < 1e-14) {
          while (min((x[-1] - x[-nsamp])/(x[nsamp] - x[1])) < 1e-06) {
              currx <- min((x[-1])[(x[-1] - x[-nsamp])/(x[nsamp] - 
                  x[1]) < 1e-06])
              k <- length(x[abs(x - currx)/(x[nsamp] - x[1]) < 
                  1e-06])
              x[abs(x - currx)/(x[nsamp] - x[1]) < 1e-06] <- currx + 
                  accuracy * (-0.5 + 1/(2 * k) + ((1:k) - 1)/k)
          }
        }
      fdist.y <- c(seq(0, 1, length = nsamp))
      }


    if (bandwidth > 0) 
        eps <- rep(bandwidth, nsamp)
    else {
        currprecision <- 0.5
        eps <- rep(0.5, nsamp)
    }
    eps[1] <- 0
    eps[nsamp] <- 0
    repeat {
        lower <- fdist.y - eps
        upper <- fdist.y + eps
        fts <- tautstring(x, fdist.y, lower, upper, upper[1], 
            lower[length(lower)], extrmean = extrema.mean)
        x.string <- fts$string
        if (sum(x.string[-1] != x.string[-(nsamp - 1)]) > 0) {
            ind1 <- min((1:(nsamp - 2))[x.string[-1] != x.string[-(nsamp - 
                1)]])
            ind2 <- max((1:(nsamp - 2))[x.string[-1] != x.string[-(nsamp - 
                1)]])
            if (x.string[ind1] > x.string[ind1 + 1]) 
                fts$nmax <- fts$nmax + 1
            if (x.string[ind2] < x.string[ind2 + 1]) 
                fts$nmax <- fts$nmax + 1
        }
        lastunif <- approx(fts$knotst, fts$knotsy, x)$y
        if (verbose) {
            par(mfrow = c(2, 1))
            if(DISCR)
              {
              plot(datax, dataemp[-1], col = "grey")
              lines(datax, x.string, col = "red")
              }
            else
              {
              hist(x, 40, prob = TRUE)
              lines((rep(x, rep(2, length(x))))[-c(1, 2 * length(x))], 
                rep(x.string, rep(2, length(x.string))), col = "red")
              }
            plot(x, upper, type = "l")
            lines(fts$knotst, fts$knotsy, col = "red")
            lines(fts$knotst, fdist.y[fts$knotsind], col = "green")
            lines(x, lower)
        }
        if (bandwidth > 0) 
            break
        if (extrema.nr > 0) {
            if (fts$nmax > extrema.nr) 
                eps <- eps + currprecision
            if (currprecision < tolerance) {
                if (fts$nmax <= extrema.nr) 
                  break
            }
            else {
                currprecision <- currprecision/2
                eps[eps > 0] <- eps[eps > 0] - currprecision
            }
        }
        else {
            diff <- cumsum(dataemp) - lastunif
            currkkuip <- kkuip(diff, maxkuipnr)$met
            kuipinds <- c(currkkuip[1], currkkuip[-1] - currkkuip[-maxkuipnr]) > 
                currbounds + 1e-08
            if (sum(kuipinds) == 0) 
                break
            eps[eps > 0] <- (currbounds[kuipinds])[1]/2
        }
        if (verbose) {
            print("Press Enter")
            readline()
        }
    }
    list(y = x.string, widthes = upper - fdist.y, nmax = fts$nmax, 
        ind = fts$knotsind, trans = lastunif)
}
"kuipdiffbounds" <-
structure(c(0.307501980728653, 0.220314152964105, 0.157498224961261, 
0.101218629121418, 0.0717806990191302, 0.0508098522102511, 0.0321599217349666, 
0.260975723191012, 0.188679698351668, 0.137057382929534, 0.0881050582317245, 
0.0633510908581634, 0.0449285678974546, 0.0285894565096712, 0.148978663288177, 
0.110940438285337, 0.080922168300817, 0.0525782956802672, 0.0374876315772114, 
0.0269628371478656, 0.0170712581613634, 0.124382620372917, 0.0952923972152614, 
0.069502435660574, 0.0452402857283996, 0.0328627738951723, 0.0234861641266666, 
0.0150209787920118, 0.0931979549885054, 0.0721989123993246, 0.0533914989167988, 
0.0353843290154652, 0.0259823805766242, 0.0185174716975161, 0.0117425052969115, 
0.080181867346804, 0.0628086817228552, 0.0476545800307207, 0.0316471006759641, 
0.0230176879538813, 0.0164825579792217, 0.010551360533595, 0.0661650592780638, 
0.0526433140545621, 0.0404992596271707, 0.0270039151375531, 0.019671300684033, 
0.0141403144918236, 0.00911192652213594, 0.0576892057496199, 
0.047291120062897, 0.0363065131491498, 0.0244300390839729, 0.017840792812333, 
0.0128944894170847, 0.0083147603774937, 0.0497274415042684, 0.0412695182684986, 
0.0319415193967662, 0.0216317776162515, 0.0159815894358178, 0.0115404374512996, 
0.0074691950384487, 0.0440930205642902, 0.0371632683985669, 0.0291764170236889, 
0.0199711130191836, 0.0147851456503776, 0.0108183394040534, 0.00692649071583678, 
0.0386375955093702, 0.0335325558345777, 0.0264992890334411, 0.0183791050868901, 
0.0136438144386473, 0.00985636827251102, 0.00638630984182286, 
0.0350372845631238, 0.0306666019917607, 0.0246776701228097, 0.0171324304025882, 
0.0127252321205716, 0.00922867868345428, 0.00599548336578164, 
0.0316423784749728, 0.0280847488931776, 0.0227942157567389, 0.0158230403649484, 
0.011929053040925, 0.00864599572624527, 0.00564937236258158, 
0.0290331314985375, 0.0259967548049147, 0.0212965842285101, 0.0149711609910706, 
0.0112445050423623, 0.0082202319675742, 0.00535064138925489, 
0.025883090604231, 0.0240572813028095, 0.0200325810700139, 0.014186955977593, 
0.0105454708236573, 0.00776751548779649, 0.00506143758819203, 
0.0232576453467689, 0.0225364463917157, 0.0188249152247203, 0.0134862219744271, 
0.0100617604237537, 0.00742707159158747, 0.00481720923365059, 
0.0206750232753329, 0.0208361558801763, 0.017644365356675, 0.0128059279425379, 
0.00958895744708681, 0.00710294102891937, 0.00459511562898805, 
0.0191363489411528, 0.0194317742601465, 0.0167332157118633, 0.0121466576268726, 
0.00920060897330519, 0.00677157289599328, 0.00443133999240373, 
0.01789137940886, 0.0183482298482047, 0.0159440097924292, 0.0116206451168149, 
0.00879431813865428, 0.0064837615919975, 0.004265070149443), .Dim = c(7, 
19))
"kuipdiffbounds.x" <-
c(50, 100, 200, 500, 1000, 2000, 5000)
"kkuip" <- function (x, k = 1)
{
    tmp <- .C("kkuip", as.double(x), as.integer(length(x)), as.integer(k),
        norm = double(k), a = integer(k), b = integer(k),PACKAGE="ftnonpar")
    list(metric = tmp$norm, a = tmp$a, b = tmp$b)
}
"rtennormal" <-
function (n)
{
    rsamp <- sample(1:10, n, rep = TRUE, prob = rep(0.1, 10))
    mus <- 0.5*(10 * rsamp - 5)
    sigmas <- 1
    rnorm(n, mus, sigmas)
}
"rclaw" <- function (n)
{
    rsamp <- sample(0:5, n, rep = TRUE, prob = c(0.1, 0.1, 0.1,
        0.1, 0.1, 0.5))
    mus <- double(n)
    sigmas <- double(n)
    mus[rsamp != 5] <- rsamp[rsamp != 5]/2 - 1
    mus[rsamp == 5] <- 0
    sigmas[rsamp != 5] <- 0.1
    sigmas[rsamp == 5] <- 1
    rnorm(n, mus, sigmas)
}
"dclaw" <- function (x)
{
    out <- 0.5 * dnorm(x)
    for (i in 0:4) out <- out + 0.1 * dnorm(x, i/2 - 1, 0.1)
    out
}

"multiwdwr" <-
function (y, thresh,firstwidth=1) 
{
    .C("multiwdwr", y = as.double(y), as.integer(length(y)), 
        as.double(thresh),as.integer(firstwidth),PACKAGE="ftnonpar")$y
}
"pmreg" <-
function (y, thr.const = 2.5, verbose = FALSE, extrema.nr = -1, 
    bandwidth = -1, sigma = -1, localsqueezing = TRUE, 
    squeezing.factor = 0.5,tolerance=0.001,extrema.mean = TRUE) 
{
 if (extrema.nr > -1) localsqueezing <- FALSE
    nsamp <- length(y)
    x<-seq(0,1,len=nsamp)
    if (sigma < 0) {
        sigma <- mad((y[-1] - y[-nsamp])/sqrt(2))
        if (verbose) 
            print(c("sigma is ", sigma))
    }
    fdist <- c(0, cumsum(y))/nsamp
    fdistx <- seq(0,1,len=nsamp+1)
    if (bandwidth < 0) 
        d <- 0.5 * (max(fdist) - min(fdist))
    else d <- bandwidth
    currprecision <- d
    eps <- rep(d, nsamp + 1)
    lower <- fdist - d
    upper <- fdist + d
    repeat {
        tstring <- tautstring(fdistx, fdist, lower, upper, 0, fdist[nsamp + 1],extrmean=extrema.mean)
        y.string <- tstring$string
        if((bandwidth<0)&&(extrema.nr<0))
          {
          residuals <- y - y.string
          residuals <- residuals - mean(residuals)
          residuals.wr <- multiwdwr(residuals, sqrt(thr.const * log(nsamp)) * sigma)
          }
        if (verbose) {
            par(mfrow = c(2, 2))
            plot(fdistx, lower, type = "l", ylim = range(c(lower, 
                upper)))
            lines(fdistx, upper)
            if (length(tstring$knotsy) > 0) {
                lines(tstring$knotst, tstring$knotsy, col = "red")
            }
            plot(x, y, col = "grey")
            lines(x, y.string, col = "red")
            if((bandwidth<0)&&(extrema.nr<0))
              plot(x, residuals.wr, type = "l", col = "green")
            print("Press Enter")
            dum <- readline()
        }
        if(bandwidth>0) break
        if(extrema.nr>0)
          {
          if(tstring$nmax>extrema.nr)
            eps<-eps+currprecision 
          if(currprecision<tolerance)
            {
            if(tstring$nmax<=extrema.nr)
            break
            }
          else
            {
            currprecision<-currprecision/2
            eps<-eps-currprecision
            }
          }
        else
          {
          ind <- (abs(residuals.wr) > 1e-10)
          ind2 <- c(FALSE, ind) | c(ind, FALSE)
          if (length(ind[ind == TRUE]) == 0) 
            break
          if (localsqueezing) 
            eps[ind2] <- eps[ind2] * squeezing.factor
          else
            eps <- eps * squeezing.factor
          }
        lower <- fdist - eps
        upper <- fdist + eps
        }
    list(y = y.string, sigma = sigma, widthes = upper - fdist, 
        nmax = tstring$nmax, knotsind = tstring$knotsind, knotsy = tstring$knotsy)
}

"tautstring" <-
function (ttt, fdist, y.low, y.up, y1 = 0.5 * (y.low[1] + y.up[1]), 
    yn = 0.5 * (y.low[length(x)] + y.up[length(x)]),extrmean=TRUE)  
{
        tmp <- .C("tautstring", as.double(fdist), as.double(ttt), 
            as.double(y.low), as.double(y.up), as.double(y1), 
            as.double(yn), as.integer(length(y.low)), tautstring = double(length(y.low) - 
                1), knotsind = integer(length(y.low)), knotst = double(length(y.low)), 
            knotsy = double(length(y.low)), nknots = integer(1), 
            nmax=integer(1),extrmean=as.integer(extrmean),PACKAGE="ftnonpar")
        list(string = tmp$tautstring, knotsind = tmp$knotsind[1:tmp$nknots], 
            knotst = tmp$knotst[1:tmp$nknots], knotsy = tmp$knotsy[1:tmp$nknots], 
            nknots = tmp$nknots,nmax=tmp$nmax)
}

"pmlogreg" <-
function (y, thr.const = 2.5, verbose = FALSE, extrema.nr = -1, bandwidth = -1, 
    localsqueezing = TRUE, squeezing.factor = 0.5, tolerance = 0.001,extrema.mean=TRUE) 
{
    if (extrema.nr > -1) localsqueezing <- FALSE
    nsamp <- length(y)
    x <- seq(0, 1, len = nsamp)
    fdist <- c(0, cumsum(y))/nsamp
    fdistx <- seq(0, 1, len = nsamp + 1)
    if (bandwidth < 0) 
        d <- 0.5 * (max(fdist) - min(fdist))
    else d <- bandwidth
    currprecision <- d
    eps <- rep(d, nsamp + 1)
    lower <- fdist - d
    upper <- fdist + d
    repeat {
        tstring <- tautstring(fdistx, fdist, lower, upper, 0, 
            fdist[nsamp + 1],extrmean=extrema.mean)
        y.string <- tstring$string
        if ((bandwidth < 0) && (extrema.nr < 0)) {
            residuals <- (y - y.string)/(sqrt(max(0.000001,y.string*(1-y.string))))
            residuals.wr <- multiwdwr(residuals, sqrt(thr.const * 
                log(nsamp)) )
        }
        if (verbose) {
            par(mfrow = c(2, 2))
            plot(fdistx, lower, type = "l", ylim = range(c(lower, 
                upper)))
            lines(fdistx, upper)
            if (length(tstring$knotsy) > 0) {
                lines(tstring$knotst, tstring$knotsy, col = "red")
            }
            plot(x, y, col = "grey")
            lines(x, y.string, col = "red")
            if ((bandwidth < 0) && (extrema.nr < 0)) 
                plot(x, residuals.wr, type = "l", col = "green")
            print("Press Enter")
            dum <- readline()
        }
        if (bandwidth > 0) 
            break
        if (extrema.nr > 0) {
            if (tstring$nmax > extrema.nr) 
                eps <- eps + currprecision
            if (currprecision < tolerance) {
                if (tstring$nmax <= extrema.nr) 
                  break
            }
            else {
                currprecision <- currprecision/2
                eps <- eps - currprecision
            }
        }
        else {
            ind <- (abs(residuals.wr) > 1e-10)
            ind2 <- c(FALSE, ind) | c(ind, FALSE)
            if (length(ind[ind == TRUE]) == 0) 
                break
            if (localsqueezing) 
                eps[ind2] <- eps[ind2] * squeezing.factor
            else eps <- eps * squeezing.factor
        }
        lower <- fdist - eps
        upper <- fdist + eps
    }
    list(y = y.string, widthes = upper - fdist, 
        nmax = tstring$nmax, knotsind = tstring$knotsind, knotsy = tstring$knotsy)
}
.First.lib <- function(lib, pkg) {
  if(version$major==0)
    stop("This version for R 1.00 or later")
  library.dynam("ftnonpar", pkg, lib)
}
"pmspec"<-function(x,pks=0,alpha=0.9,sqzf=0.9,mult=0,lcl=0,ln=0,fig=0,pow=10^-2)  
{
        m<-length(x)
	kk<-1	
	while (2**kk <m ) {
			 kk<-kk+1
			 }
	n<-2**kk
	edf<-double(n)
	edf[1:n]<-0		
	edf[1:m]<-x-mean(x)
	edf<-fft(edf)
	edf<-Re(edf)**2+Im(edf)**2
	edf<-edf/(4*pi*m)	
	n<-n/2
	edf<-edf[1:n]
	if(n<=256){
		mult<-1
		}
	edf<-edf+pow*mean(edf)
	kk<-kk-1
#
#
	sqzf<-min(0.95,sqzf)	
#	calculate cutoff value
#	beta2<-1-(1-alpha)/(2*n)
#	beta1<-1-beta2
	if (mult > 0){
		kk<-n-1
		}
	qxtrm<-double(2*kk+2)
#
                ### fortran subroutine###
#	
    	tmp <- .Fortran(
		"npspcdn",
		as.double(edf),
		double(n+1),
		double(n+1),
		double(n+1),
		double(n),
		double(n),
		double(n+1),
		as.double(qxtrm),
		integer(n+1),
		integer(n+1),
		integer(n+1),
		as.integer(n),
		integer(1),
		as.double(sqzf),
		as.integer(pks),
		as.integer(kk),
		as.integer(mult),
		as.integer(lcl),
		as.double(alpha),
                PACKAGE="ftnonpar"
		)
#	
              ###spectral density###
#
	df<- tmp[[5]]
	edf<-edf
	sy<-tmp[[2]]
	ll<-tmp[[3]]
	uu<-tmp[[4]]
	nkn<-tmp[[12]]
	str<-double(2*nkn)
	dim(str)<-c(nkn,2)
	dm<-tmp[[9]]
	str[,1]<-dm[1:nkn]
	dm<-tmp[[7]]
	str[,2]<-dm[1:nkn]
	pks<-tmp[15]
#
              ###number of peaks##
#		
#              ###plot scale###
#
	if(fig==0){	
                xx<-0:(n-1)
		xx<-pi*xx/n
	if(ln==0) {
		  low<- min(log(edf[2:(n-1)]),log(df[2:(n-1)]))		
	          upp<- max(log(edf[2:(n-1)]),log(df[2:(n-1)]))
		  rng<-upp-low
		  low<-low-0.02*rng
		  upp<-upp+0.02*rng
	          plot(xx[2:(n-1)],log(edf[2:(n-1)]),ylim=range(low,upp),
	              xlab="Frequency ",ylab="Log(spectral density)",
			col=2) 
	          lines(xx[2:(n-1)],log(df[2:(n-1)]))
	          }
	else      {
	          upp<- 1.02*max(edf,df)
	          plot(xx,df,type="l",ylim=range(0,upp),
	              xlab="Frequency ",ylab="Spectral density") 	
	          points(xx,edf,col=2)
	          }
	}
#
	list(edf=edf,df=df,pks=pks,ll=ll,uu=uu,str=str,qxt=qxtrm)
	
#  
#
#
                       ###INPUTS###
#
# x	     = data
#
# alpha      = level for scaled observations
#
# sqzf	     = squeeze factor for tube
#
	               ###OUTPUTS###
#
# scl	     = scale
#
# pks	     = number of peaks
#							  	
}


