program protein;  { Computes and saves polypeptide extensions for many
log-spaced forces using power-series WLC from Bouchiat et al,
Biophysical J, 1999. }

{$N+}
uses crt,graph;
const
  kB = 1.38064e-23; {Boltzmann const in MKS}
  P1 = 6.5e-10; {persistence length of polypeptide at 25 deg.C}
var snert,tension,f,f2,x,lastX,dX,L,P,T,temp:real;
    ch:char;
    i,j,n:integer;
    s:array[1..7] of real;
    outFile:text;
    message:string;
begin
clrscr;  writeln;
writeln('What to call the output file? Extension .txt will add automatically');
readln(message);
assign(outFile,message + '.txt');
rewrite(outFile);
writeln('what temperature deg.C ?');
readln(temp);   T := temp + 273.15;
write(outFile,'generated by Protein.pas    T= ', temp:1:1 ,' deg.C');
P := P1 * (273.15 + 25) / T;    {P = bending rigidity / kT}
writeln(outFile,'  P = ',P*1e9:1:3,' nm');
writeln(outFile);
writeln(outFile,'force(pN)             fractional extension');
tension := 4e-12; {Newtons}
j:=1;

repeat
snert:=kB*T/(4*P);
L:=(1);
x:=0.5;  n:=1;
repeat
      s[1]:=x/L;
      for i:=2 to 7 do s[i]:=s[1]*s[i-1];
      f := snert*(1/sqr(1-(x/L))+4*(x/L)-1-2.0656912*s[2]-10.949672*s[3]
      +64.29988*s[4]-155.50428*s[5]+157.99776*s[6]-56.70872*s[7]);
      f2 :=(2*snert/L)*(1/((1-(x/L))*(1-(x/L))*(1-(x/L)))+2);
      if (tension-f)/(n*f2) < (L-x) then x:=x+(tension-f)/(n*f2) else inc(n);
until abs(f-tension)<0.01*tension;
write(outFile,tension*1e12:0:3);
writeln(outFile,x:30:3);     {fractional extension}
tension:=tension*1.005;
until (tension>4e-11);
close(outFile);
end.