(*Computation of values for canonical sec str  Nov 23, 2008 *)
 Clear[R,S,P,Q,F,K1,K2];



eqn = {S == z + S z +R z^2 + R S  z^2,R ==  z^3 + S z^3 +  R S z^4 }
sol = Solve[eqn,S,R]

sol1 =Last[ Extract[sol,{1}]]

Series[4 (z^5) (-1-(z^4))+(((-1+z+(z^5)))^2),{z,0,40}]

F =1-2 z+(z^2)-6 (z^5)+2 (z^6)-4 (z^9)+(z^10)

solset=NSolve[F==0, z]

rho =0.5309527613419542`

(* dominant singularity *)

(* Now rho is dominant singularity.Want to factor (1-z/rho) from P.*)

growth = 1/rho
Q=Series[F/(1-z/rho),{z,0,10}]

Print[F]
Expand[Series[Q,{z,0,10}]*(1-z/rho)]

(*So F=Q(1-z/rho) *)

(* Q0 is approximation to Q, which is F with (1-z/rho) factored out*)
Q0 =1-0.11659327757792548` z+0.7804074372205121` (z^2)+1.469824613489295` (z^3)+2.768277557827167` (z^4)-0.7862074380581525` (z^5)+0.5192516259430491` (z^6)+0.9779620029297007` (z^7)+1.8419002105912057` (z^8)-0.5309527613419505` (z^9)
K0=Q0/. z->rho

K1 = -Sqrt[K0]/(2 (z^4))/.z->rho

K2 =K1/Gamma[-1/2]

(* asymptotic function *)

K[n_] := K2*n^(-3/2)*(1/rho)^n
For[i=1,i≤50,i++,Print[i,"  ",K[i]]]

(* Compute average number of base pairs in sat str, where theta is 0*)
s = 1/(1-u z^2 )
Series[D[s,u]/.u->1,{z,0,30}]
Series[s,{z,0,50}]

(********************************************
Canonical structures, comparison with Hofacker et al;
Here we get the same dominant singularity as in Table 1 on page 15 of Hofacker et al. where value of L is 2,minimum stack size, and value of m is 1,minimum number of unpaired bases in hairpin loop.
************************************************)


-S+((1-z-(z^5)-(Sqrt[4 (z^5) (-1-(z^4))+(((-1+z+(z^5)))^2)]))/(2 (z^4)))=0
S*2z^4 -1+z+z^5+(Sqrt[4 (z^5) (-1-(z^4))+(((-1+z+(z^5)))^2)]) =0
F=Expand[(S*2z^4 -1+z+z^5+(Sqrt[4 (z^5) (-1-(z^4))+(((-1+z+(z^5)))^2)]))^2]
Fx = D[F,z]=-4 S (z^3) (-S+z+S z)+((z^2)+S (z^2)) (3 (z^2)+3 S (z^2))+(2 z+2 S z) ((z^3)+S (z^3))-(1+S) (-1+S (z^4));
Fyy = D[D[F,S],S]=-2 (-1+z) (z^4)+2 (z^5);
alpha = 0.5081
beta = (1/alpha^2)*Sqrt[1-alpha^2+alpha^4]
FxEval = Fx/. {z->alpha,S->beta}
FyyEval = Fyy /. {z->alpha,S->beta}
gAlpha = -Sqrt[(alpha * FxEval)/(2*Pi*FyyEval)]
growthRate = 1/alpha
constant = -gAlpha/(2Sqrt[Pi])


F=Resultant[-S+z+S z + R z^2 + S R z^2,-R + z^3 + R z^2 + S R z^4 + S z^3,R];
F=((z^2)+S (z^2)) ((z^3)+S (z^3))-(-S+z+S z) (-1+(z^2)+S (z^4));
Fx = D[F,z]=((z^2)+S (z^2)) (3 (z^2)+3 S (z^2))+(2 z+2 S z) ((z^3)+S (z^3))-(-S+z+S z) (2 z+4 S (z^3))-(1+S) (-1+(z^2)+S (z^4));
Fyy = D[D[F,S],S]=-2 (-1+z) (z^4)+2 (z^5);
alpha = 0.5081
beta = (1/alpha^2)*Sqrt[1-alpha^2+alpha^4]
FxEval = Fx/. {z->alpha,S->beta}
FyyEval = Fyy /. {z->alpha,S->beta}
gAlpha = -Sqrt[2alpha * FxEval/FyyEval]
growthRate = 1/alpha
constant = -gAlpha/(2Sqrt[Pi])

Eliminate[{S==z+S z + R z^2 + S R z^2,R== z^3 + R z^2 + S R z^4 + S z^3},R]
F=(S^2) (z^4)+S (-1+z+(z^2)-(z^3)+(z^5))-z (-1+(z^2)-(z^4))
Fx = D[F,z]
Fyy = D[D[F,S],S]
alpha = 0.5081
beta = (1/alpha^2)*Sqrt[1-alpha^2+alpha^4]
FxEval = Fx/. {z->alpha,S->beta}
FyyEval = Fyy /. {z->alpha,S->beta}
gAlpha = -Sqrt[2alpha * FxEval/FyyEval]
growthRate = 1/alpha
constant = -gAlpha/(2Sqrt[Pi])

(**************************************************************
By taking F defined by solution of equations or resultant or from elimination of variables I get in both cases that constant is 2.7437038542949987` rather than the real value 2.1614.
**************************************************************)