"corrlen"<-
function(sl = 0.75, sh = 1.25, d = 0.01, sigma = cursigma)
{
	slist <- seq(sl, sh, d)
	l <- numeric()
	for(s in slist) {
		if(s < 1)
			l <- c(l, 1/(1 - s))
		if(s == 1)
			l <- c(l, s/sqrt(sigma))
		if(s > 1)
			l <- c(l, (s - 1)/sigma + 1/(s - 1))
	}
	cbind(slist, l)
}
"dlnqdz"<-
function(s = seq(0.90000000000000002, 1.1000000000000001, 0.02), o = cursigma, ds
	 = 1, do = 0, n = 1000)
{
# computation of derivatives of partition function, Q;
#
# 1-d Ising model of helix-coil transition (Zimm and Bragg, 1959);
#
# for mu = 1, where mu is the number of coil segments that must precede
#	the helix segment at the start of a helix
#
# input parameters:
#	s = single value or sequence of values of s, the factor
#		contributed by a helix segment to the statistical weight 
#		of a chain configuration
#	o = value of sigma, the additional factor contributed to the
#		statistical weight by the difficulty of initiating
#		a helical region
#	ds = 1 if the derivative returned is dlnQ/ds; ds = 0 otherwise
#	do = 1 if the derivative returned is dlnQ/dsigma; do = 0 otherwise
#	n = length of (number of segments in) the chain
#
# value returned:
#	vector:
#		cbind(s, dlnq, sdlnq, odlnq, l0, l1, theta, dthdlns)
#
#		elements:
#			s
#			dlnq	= dlnQ/dz, where dz == ds or do
#			sdlnq	= s * dlnQ/ds
#				= <k>	= number of helix segments per chain
#			odlnq	= sigma * dlnQ/dsigma
#				= nu	= number of helices per chain
#			l0	= larger eigenvalue
#			l1	= smaller eigenvalue
#			theta	= (s/(n-3))*dlnQ/dz
#				= fraction helix
#			dthdlns	= dtheta/dlns
#				= sharpness of transition
#
#
	D <- sqrt((1 - s) * (1 - s) + 4 * o * s)
	C <- ( - ds * (1 - s) + 2 * o * ds + 2 * s * do)/D
	l0 <- (1 + s + D)/2
	l1 <- (1 + s - D)/2
	dl0 <- (ds + C)/2
	dl1 <- (ds - C)/2
	N1 <- (l1 - 1) + (1 - l0) * (l1/l0)^(n - 2)
	Denom <- l1 - l0
	dN1 <- (n - 2) * dl0 * (l1 - 1) + l0 * dl1 + ((n - 2) * dl1 * (1 - l0) - 
		l1 * dl0) * (l1/l0)^(n - 3)
	dDenom <- (dl1 - dl0)
	dlnq <- dN1/(l0 * N1) - dDenom/Denom
	dthdlns <- rep(1, length(s))
	if(ds == 1) {
		ddl0 <- (1 - C * C)/(2 * D)
		ddl1 <-  - ddl0
		A <- (n - 2) * (ddl0 * (l1 - 1) + dl0 * dl1) + dl0 * dl1 + l0 * 
			ddl1
		B <- ( - l1 * dl0)/(l0 * l0) + dl1/l0
		A <- A + ((n - 2) * dl1 * (1 - l0) - l1 * dl0) * (l1/l0)^(n - 4
			) * (n - 3) * B
		A <- A + (l1/l0)^(n - 3) * ((n - 2) * (ddl1 * (1 - l0) - dl1 * 
			dl0) - dl1 * dl0 - l1 * ddl0)
		A <- (s * A)/(l0 * N1)
		X <- dl0 * ((l1 - 1) + (l1/l0)^(n - 2) * (1 - l0))
		X <- X + l0 * (dl1 - (l1/l0)^(n - 2) * dl0)
		X <- X + l0 * (n - 2) * (1 - l0) * (l1/l0)^(n - 3) * B
		X <- ( - s * X * dN1)/(l0 * N1 * l0 * N1)
		Y <- (dl1 - dl0)^2/(l1 - l0)^2 - (ddl1 - ddl0)/(l1 - l0)
		Y <- s * Y
		dthdlns <- s/(n - 3) * (A + X + Y + dlnq)
	}
	sdlnq <- s * dlnq
	odlnq <- o * dlnq
	theta <- sdlnq/(n - 3)	

#	del <- diff(theta)/diff(s)
#	del <- c(del[1], del)
#	f <- s/(n - 3)
#	cbind(s, dlnq, sdlnq, odlnq, l0, l1, theta, del, dthdlns, A = f * A, X
#		 = f * X, Y = f * Y, q = f * dlnq, ddl0, ddl1)

	cbind(s, dlnq, sdlnq, odlnq, l0, l1, theta, dthdlns)
}
"fnh"<-
function(s = seq(0.75, 1.25, 0.01), sigma = cursigma, n = 1000)
{
	kav <- dlnqdz(s, sigma, 1, 0, n)[, "sdlnq"]
	nu <- dlnqdz(s, sigma, 0, 1, n)[, "odlnq"]
	lav <- kav/nu
	theta <- kav/(n - 3)
	cbind(s, lav, theta, kav, nu)
}
"fnhelix"<-
function(sl = 0.75, sh = 1.25, d = 0.01, sigma = cursigma)
{
	slist <- c(seq(sl, 0.97999999999999998, d), 1, seq(1.02, sh, d))
	theta <- numeric(0)
	for(s in slist) {
		if(s < 0.98999999999999999)
			theta <- c(theta, (sigma * s)/((1 - s + sigma * s) * (1 -
				s)))
		else if(s > 1.01)
			theta <- c(theta, (s - 1 - sigma/(s - 1))/(s - 1 - 
				sigma))
		else if(s == 1)
			theta <- c(theta, (s/2)/((1 + s)/2 + sqrt(sigma)))
	}
	cbind(slist, theta)
}
"fnmat"<-
function(sigma = cursigma)
{
	slast <- 0.10000000000000001
	s <- slast
	slist <- for(i in 1:40) {
		slast <- slast * 10^(1/20)
		s <- c(s, slast)
	}
	s <- slist
	nlast <- 4
	n <- nlast
	nlist <- for(i in 1:50) {
		nlast <- nlast * 10^(1/10)
		n <- c(n, nlast)
	}
	for(n in nlist) {
		kav <- dlnqdz(s, sigma, 1, 0, n)[, "sdlnq"]
		nu <- dlnqdz(s, sigma, 0, 1, n)[, "odlnq"]
		theta <- kav/(n - 3)
		if(n == nlist[1])
			nblock <- cbind(s, n, theta, nu)
		else nblock <- rbind(nblock, cbind(s, n, theta, nu))
	}
	list(s = s, n = nlist, data = nblock)
}
"newmat"<-
function(sigma = cursigma)
{
	slast <- 0.10000000000000001
	s <- slast
	slist <- for(i in 1:40) {
		slast <- slast * 10^(1/20)
		s <- c(s, slast)
	}
	s <- slist
	nlast <- 4
	n <- nlast
	nlist <- for(i in 1:50) {
		nlast <- nlast * 10^(1/10)
		n <- c(n, nlast)
	}
	for(n in nlist) {
		temp <- dlnqdz(s, sigma, 1, 0, n)
		kav <- temp[, "sdlnq"]
		dtheta <- temp[, "dthdlns"]
		theta <- temp[, "theta"]
		nu <- dlnqdz(s, sigma, 0, 1, n)[, "odlnq"]
		if(n == nlist[1])
			nblock <- cbind(s, n, theta, nu, dtheta)
		else nblock <- rbind(nblock, cbind(s, n, theta, nu, dtheta))
	}
	th <- matrix(nblock[, "theta"], nrow = 41, ncol = 51, byrow = F)
	dth <- matrix(nblock[, "dtheta"], nrow = 41, ncol = 51, byrow = F)
	numat <- matrix(nblock[, "nu"], nrow = 41, ncol = 51, byrow = F)
	list(s = s, n = nlist, data = nblock, th = th, dth = dth, nu = numat)
}