initialise:=
Block[{},
hh[r_,s_] := (r^2-1)/(4t) + (s^2-1)t/4 - (r s-1)/2
;
Clear[head,tail,a,b,c,cj,xj,tot]
;
Clear[Amj];
Amj[ a___ ] := Apply[ mode, Join[ a ]]
;
Clear[App];
App[{a_,n_},exp___] := Block[{c},
  exp /. mode[c___] -> mode[{a,n},c] ]
  ;
Clear[Ppa];
Ppa[{a_,n_},exp___] := Block[{c},
  exp /. mode[c___] -> mode[c,{a,n}] ]
  ;
Clear[mode];
setvac:= Block[{c,a},
   mode/: mode[c___, {L,a_}] := 0 /; a > -2 ;   ]
   ;
sethwt:= Block[{c,a},
   mode/: mode[c___, {L,0}]  := h mode[c] ;
   mode/: mode[c___, {L,a_}] := 0  /; a > 0 ; ]
;
mode/: mode[c___,{L,a1_},{L,a2_},b___] :=  
  Block[{},  mode[c,{L,a1},{L,a2},b] =  
    Expand[ mode[ c, {L,a2}, {L,a1}, b] + 
            CommLL[ a1 , a2 , {c}, {b} ]  ]
	] /; a1 > a2 
	; 
Clear[CommLL]
;
CommLL[ a1_,a2_,a3_List,a4_List] :=  Block[ {},
    If[a1 + a2 ==0,  cc*(a1 -1)*a1 *(a1 +1)/12,0] * Amj[ a3, a4 ]  +
    (a1 - a2) * Amj[ a3, {{L,a1 + a2 }}, a4 ] 
    ]
    ;
Clear[Changefield];
Changefield[L]=l;
mkmode :=
 Block[{},
  Clear[mode];
  mode[a___] := 
   Block[{b,ans},
   ans = {a};
   fn1 = Changefield[#[[1]]][- #[[2]] ]& ;
   ans = Map[  fn1, List @@ ans ];
  ans = NonCommutativeMultiply @@ #& /@  {ans} ; 
  ans =  # ** ket[0,0]& @@ #& /@ {ans} ;
   Return[ans];
	]
]
;
Clear[basis];
basis[0] = { mode [] } ;
basis[j_] := basis[j] = addbasis[j,mode[]];
Clear[addbasis];
addbasis[ target_, vec_mode] :=
Block[{k, ans, sofar, remainder, lvec, k1, k2}, 
    If[ TrueQ[vec == mode[] ], 
	  k1 = X;
	  k2 = -target;
	  sofar = 0;
	  remainder = target;
	, lvec = Last[vec];
	  k1 = lvec[[1]];
	  k2 = lvec[[2]];
	  sofar = - (Plus @@ ( #[[2]]& /@ vec ) ) ;
          remainder = target - sofar;
	, Print[ "???" ]
	];
    If[ k1 == X , 
	  ans = {
	      addL[ -remainder,target,vec ]
	      }
	  ];
    If[ k1 == L , 
	  ans = {
	      addL[  Max[k2,-remainder],target,vec ]
	      }
	  ];
    If[ remainder == 0, ans = vec ];
      ans = Flatten[ans]; 
      Return[ans]
      ];
Clear[addL]  ;
 ;   
addL[start_,target_,vec_mode]  :=
    addbasis[target,#]& /@
    Block[{k},
    Table[
	Join[ vec, mode[{L,k}] ],
	{k, -Floor[-start], -1}
	]];
Clear[m1,m2];
m1[level_]  := m1[ level ] =
   Block[ {exp},
       b1 = basis[level ];
       b2 = basis[level - 1 ];
       exp = App[{L,1}, b1 ];
       l2 = Length[b2] ;
      Table[ D[exp, b2[[j]] ],{j,1,l2} ]
      ];
    m2[level_]  := m2[ level ] =
   Block[ {exp},
       b1 = basis[level ];
       b2 = basis[level - 2 ];
       exp = App[{L,2}, b1 ];
       l2 = Length[b2] ;
      Table[ D[exp, b2[[j]] ],{j,1,l2} ]
      ];
findnv0[j_ ]:= Join[ m1[j],m2[j]];    
Clear[findnv];  
findnv[j_ ]:= findnv[j] = NullSpace[ findnv0[j]] ;
;
;
Clear[D1,D0,DL];
D1[ f_ ]:= - z D[f,z]           ;
D0[ f_ ]:= (z-1)D[f,z]          ;
DL[n_, f_ ] := -( 
   (1-z)^(n+1) D1[f] + 
   (0-z)^(n+1) D0[f] )          ;
Clear[dmode];
;
dmode[ {L,n_} , rest___] := 
dmode[ {L,0} , rest] + 
( (z - z^(n+1)) * D[ Evaluate[ dmode[rest]],z]
 + ( (1 - (n+1)*z^n) * h1  * dmode[rest]
 +  ( -n * h2) * dmode[rest] )
 + (-h3 * dmode[rest])  ) /; n < 0
 ; 
dmode[ {L,0}, rest___] := (Evaluate[ mode[{L,0},rest] ])/. mode -> dmode  
;
dmode[ {L,0} ] = h * y[z];
 ; 
dmode[ ] = y[z] ;
Clear[mkdiff];
mkdiff[state_] := Expand[state /. mode -> dmode ]
;
nv[level_] := (findnv[level] . basis[level])[[1]] ;

diffeq[level_] := mkdiff[ nv[level]] ;
];
(* Recursion formula for conformal block by I Runkel from
   hep-th/0107118
   with corrections from
[1] Al.B. Zamolodchikov, Conformal Symmetry In Two-Dimensions: An Explicit Recurrence Formula For The Conformal Partial Wave Amplitude, Commun. Math. Phys. 96 (1984) 419–422.
[2] Al.B. Zamolodchikov, Conformal Symmetry in Two-Dimensional Space: Recursion Representation of Conformal Block, Theor. Math. Phys. 73 (1987) 1088–1093.
 *)

     ft[d_,m_,n_] := ( 2d + m*n - 1 +
         Sqrt[ 4d(d+m*n-1) + (m-n)^2 ] ) / (n^2-1);

     fc[t_] := 13 - 6(t+1/t);

     fl[t_,st_,p_,q_] := 1/(2*st) * (p - q*t);

     fld[d_,t_] := Sqrt[ d + (1-t)^2/(4*t) ];

     fA[d_,m_,n_] := Block[{res,t,st,p,q,k,l},
          t = ft[d,m,n]; st = Sqrt[t];
          res = 12*(t-1/t) / ( (m^2-1)/t - (n^2-1)*t ) *
          Product[ If[ (k==0 && l==0) || (k==m && l==n), 1, 
               1/(2*fl[t,st,k,l]) ],
               {k,-m+1,m}, {l,-n+1,n} ];
          (*Factor[Simplify[Factor[res]]]*)
          res
          ];

     fP[dv_,d_,m_,n_] := Block[{l1,l2,l3,l4,res,t,st,p,q,k,l},
          t = ft[d,m,n]; st = Sqrt[t];
          l1 = fld[dv[[1]],t]; l2 = fld[dv[[2]],t]; 
          l3 = fld[dv[[3]],t]; l4 = fld[dv[[4]],t]; 
          res = Product[
               ( l2 + l1 - fl[t,st,p,q] ) *
               ( l2 - l1 - fl[t,st,p,q] ) *
               ( l3 + l4 - fl[t,st,p,q] ) * 
               ( l3 - l4 - fl[t,st,p,q] ) ,
               {p,-m+1,m-1,2}, {q,-n+1,n-1,2} ];
          (*Factor[Simplify[Factor[res]]]*)
          res
          ];

     FC[c_,dv_,d_,x_,lvl_] := fF[c,dv,d,x,lvl];
     fF[c_,dv_,d_,x_,lvl_] := Block[ {m,n},
          If[lvl<=0,1,

          Hypergeometric2F1[d+dv[[2]]-dv[[1]],d+dv[[3]]-dv[[4]],2d, x] + 

          Sum[ Sum[
               x^(m*n) / ( c - fc[ ft[d,m,n] ] ) *
               fA[d,m,n] * fP[dv,d,m,n] *
               fF[ fc[ft[d,m,n]], dv, d + m*n, x, lvl-m*n ] ,
               {m,1,Floor[lvl/n]}] , {n,2,lvl} 
               ]
          ] ];
(* Recursion formula for conformal block by I Runkel (hep-th/0107118) with changes from from
[1] Al.B. Zamolodchikov, Conformal Symmetry In Two-Dimensions: An Explicit Recurrence Formula For The Conformal Partial Wave Amplitude, Commun. Math. Phys. 96 (1984) 419–422.
[2] Al.B. Zamolodchikov, Conformal Symmetry in Two-Dimensional Space: Recursion Representation of Conformal Block, Theor. Math. Phys. 73 (1987) 1088–1093.
[3] A.B.Zamolodchikov & Al.B.Zamolodchikov, Structure Constants and Conformal Bootstrap in Liouville Field Theory, hep-th/9506136, section 7

 zaF[c,{h1,h2,h3,h4},h,x,level]

level = level to which to descend in the recursion
 *)

     zaq[x_] := Exp[ Pi*I*zatau[x] ];

     zatau[x_] := I * zaK[1-x] / zaK[x];

     zaK[x_] := Block[ {t}, 1/2 * 
          NIntegrate[ ( t(1-t)(1-x*t) )^(-1/2) , {t,0,1} ] ];

     zaK[x_] := EllipticK[x];

     (* alpha_plus and alpha_minus *)
     Clear[zaap,zaam];
     zaap[c_] := zaap[c] = 1/Sqrt[24] * ( Sqrt[1-c] + Sqrt[25-c] );
     zaam[c_] := zaam[c] = 1/Sqrt[24] * ( Sqrt[1-c] - Sqrt[25-c] );

     (* relation lambda <-> delta *)
     zalami[c_,dv_,n_] := Sqrt[dv[[n]]+(1-c)/24];

     (* lambda_pq *)
     zalam[c_,p_,q_] := zaap[c]*p + zaam[c]*q;

     (* R_mn *)
     zaR[c_,dv_,m_,n_] := Block[ {l1,l2,l3,l4,p,q,k,l},
          l1 = zalami[c,dv,1]; l2 = zalami[c,dv,2];
          l3 = zalami[c,dv,3]; l4 = zalami[c,dv,4];

          -1/2 * 
          Product[
               ( l2 + l1 - zalam[c,p,q]/2 ) *
               ( l2 - l1 - zalam[c,p,q]/2 ) *
               ( l3 + l4 - zalam[c,p,q]/2 ) * 
               ( l3 - l4 - zalam[c,p,q]/2 ) ,
               {p,-m+1,m-1,2}, {q,-n+1,n-1,2} ] *
          Product[ 
               If[ (k==0 && l==0) || (k==m && l==n), 
                    1, 
                    1/ zalam[c,k,l] 
               ],
               {k,-m+1,m}, {l,-n+1,n} ]
     ];

     (* Delta_mn *)
     zaDel[c_,m_,n_] := (c-1)/24 + (zaap[c]*m + zaam[c]*n)^2 / 4;

     (* recusive function H *)
     zaH[c_,dv_,del_,q_,lvl_] :=  Block[ {m,n},
          If[lvl<=0,1,
               1 +           
               Sum[ Sum[
                    (16*q)^(m*n) / ( del - zaDel[c,m,n] ) *
                    zaR[c,dv,m,n] * zaH[ c, dv, zaDel[c,m,n]+m*n, q , lvl-m*n] ,
               {m,1,Floor[lvl/n]}] , {n,1,lvl} 
               ]

          ]
     ];

     (* theta3 *)
     theta3[q_] := EllipticTheta[3,0,q];

     (* conformal block F *)

     zaF[c_,dv_,d_,x_,lvl_] := Block[ {q},
          q = zaq[x];
          (16*q)^( d-(c-1)/24 ) *
          x^( (c-1)/24 - dv[[3]] - dv[[4]] ) *
          (1-x)^( (c-1)/24 - dv[[2]] - dv[[3]] ) *
          theta3[q]^( (c-1)/2 - 4*(dv[[1]]+dv[[2]]+dv[[3]]+dv[[4]]) ) *
				      zaH[c,dv,d,q,lvl]];

FH[c_,dv_,d_,x_,lvl_] := Block[{},x^( dv[[3]] + dv[[4]] -d ) zaF[c,dv,d,x,lvl]];

There is a hidden initialisation cells with definitions above this line.
It defines the commands
 “initialise” for the null-vector and differential equation section.
 “FC“ and “FH” for the Zamolodchikov recursion formulae (courtesy of Ingo Runkel)

Virasoro singular vectors and differential equations

Ising free fermion

h = 1/2 ;
cc = 1/2 ;

initialise;

sethwt

findnv[2]

tosolve = diffeq[2] /. {h1 ->h, h2 -> h,h3 -> h}

DSolve[tosolve==0,y[z],z]

Simplify[tosolve /. y -> Function[x,(1-x+x^2)/((1-x)x)]]

Ising spin field

h = 1/16; cc = 1/2 ;

initialise;

sethwt

findnv[2]

tosolve = diffeq[2] /. {h1 -> 1/16, h2 -> 1/16,h3 -> 1/16}

DSolve[tosolve==0,y[z],z]

Simplify[tosolve /. y -> Function[x,(x(1-x))^(-1/8)Sqrt[1 + Sqrt[1-x]]]]
Simplify[tosolve /. y -> Function[x,(x(1-x))^(-1/8)Sqrt[1 - Sqrt[1-x]]]]

Lee Yang

SetDirectory["/Users/gerard/projects/blocks/"]

h = -1/5; cc = -22/5 ;

<< "VirasoroDifferentialEquations.m";

sethwt

findnv[2]

tosolve = diffeq[2] /. {h1 ->h, h2 -> h,h3 -> h}

DSolve[tosolve==0,y[z],z]

FullSimplify[%]

FullSimplify[tosolve /. y -> Function[x, x^(2/5)(1-x)^(2/5)Hypergeometric2F1[3/5,4/5,6/5,x]]
]

FullSimplify[tosolve /. y -> Function[x, x^(1/5)(1-x)^(2/5)Hypergeometric2F1[2/5,3/5,4/5,x]]
]

Generic level 2

SetDirectory["/Users/gerard/projects/blocks/"]

h = 3 t/4 - 1/2; cc = 13 - 6 t - 6/t;

initialise;

sethwt

findnv[2]

tosolve = Simplify[diffeq[2] /. {h1 -> hh[1,α], h2 -> hh[1,β], h3-> hh[1,γ]} ]

DSolve[tosolve==0,y[z],z]

FullSimplify[PowerExpand[FullSimplify[DSolve[tosolve==0,y[z],z]]] /. C[2] -> 0]

FullSimplify[PowerExpand[FullSimplify[DSolve[tosolve==0,y[z],z]]] /. C[1] -> 0]

Generic level 3

SetDirectory["/Users/gerard/projects/blocks/"]

h = hh[1,3] 
cc = 13 - 6 t - 6/t;

<< "VirasoroDifferentialEquations.m";

sethwt

findnv[2]

findnv[3]

tosolve = diffeq[3] /. {h1 -> h, h2 -> h,h3 -> h}

DSolve[tosolve==0,y[z],z]

DSolve[(tosolve /. t -> 3/4 /. h -> 1/2 )== 0,y[z],z]

FullSimplify[tosolve /. y -> Function[x,(1-x+x^2)/(x(1-x))]]
% /. t -> 3/4

Zamolodchikov recursion formulae

The h-recursion formula is called by
FH[c,{h1,h2,h3,h4},h,z,level]
where level is the depth required. 
It does not include the prefactor (z^((h-h3-h4)))

Example: The free fermion block

IFFB2 = FH[1/2,{1/2,1/2,1/2,1/2},0,z,2]

To find its series expansion. This is only reliable up to a certain order.

Series[ IFFB2,{z,0,4}]

The c-recursion formula is called by
FC[c,{h1,h2,h3,h4},h,z,level]
where level is the depth required. 
It does not include the prefactor (z^((h-h3-h4)))
It is not well-defined on minimal model values where there are poles in c

Example: The free fermion block

IFFB2 = FC[1/2,{1/2,1/2,1/2,1/2},0,z,3]

They agree on numerical values and converge well, and agree

FH[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,1]
FH[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,2]
FH[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,3]
FH[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,10]

FC[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,1]
FC[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,2]
FC[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,3]
FC[100.1,{1.2,1.3,1.4,1.5},1.6,0.7,10]

Plots from the lecture, and a few more

Clear[h]
h[r_,s_,t_] := (r^2-1) /(4 t) + (s^2-1)t/(4) - (r s -1)/2

c[t_] := 13 - 6 t - 6/t

Clear[eqn]
eqn[r_,s_] := eqn[r,s] = Simplify[Expand[(h-h[r,s,t])(h - h[s,r,t])] /. t->(1/12) (13-CC-(Sqrt[25-26 CC+(CC^2)]))]

oo[r_,s_] := ContourPlot[ eqn[r,s]==0,{h,-15,15},{CC,-40,60}]

oo2[r_,s_] := Show[ParametricPlot[ {h[r,s,t],c[t]},{t,-10,10},
PlotRange -> {{-15,15},{-40,60}},AspectRatio->1],
ParametricPlot[ {h[r,s,Exp[I t]],c[Exp[I t]]},{t,0,2 Pi},
PlotRange -> {{-15,15},{-40,60}},AspectRatio->1]]

kacplot[n_] := kacplot[n] = Show[ Table[ oo2[r,s],{r,1,n},{s,1,n/r}]]

kacplot2[n_,hr_,cr_] := Show[kacplot[n],PlotRange -> {hr,cr}]

commands  = {Plot[100,
{h,1,2},PlotRange -> {{-15,15},{-40,60}},AspectRatio->1],
kacplot[1],kacplot[2],kacplot[3],kacplot[4],kacplot[10],kacplot[20],kacplot[100],
kacplot2[100,{-1,2},{0,1}],
Show[kacplot2[100,{-1,2},{0,1}],Plot[{(1/2),(7/10),(4/5),(6/7),(25/28),(11/12),(14/15),(52/55)},{h,-1,2},PlotStyle->Red],kacplot2[100,{-1,2},{0,1}]],
kacplot2[100,{-1,2},{-20,1.1}],
Show[kacplot2[100,{-1,2},{-20,1.1}],Plot[{-(22/5),-(68/7),-(46/3),-(232/11),-(350/13),-(164/5),-(658/17),-(848/19)},{h,-1,2},PlotStyle->Red]]
,
kacplot2[100,{-15,15},{0,60}]
};

Manipulate[Show[commands[[n]]],{n,1,13,1},ControlPlacement->Right]

hctab = Table[{h[r,s,t],c[t]},{r,1,20},{s,1,r}];

lines = Table[{(1/2)-((r^2)/2)-(Cosa/2)+(1/2) (r^2) Cosa,13-12 Cosa},{r,1,50}];

lineplot = ParametricPlot[ lines,{Cosa,-20,20}];

allcplot = Show[lineplot,ParametricPlot[ hctab,{t,-10,10},
PlotRange -> {{-15,15},{-40,60}},AspectRatio->1],PlotRange -> {{-15,15},{-40,60}},AspectRatio->1];

plot1 = Show[allcplot, Plot[{1,25,1 + 24 h},{h,-15,15},PlotStyle->Dashed]]

plot2 = Show[plot1,PlotRange -> {{-1,15},{0,1}},AspectRatio -> 1]

posposcplot = ParametricPlot[ hctab,{t,0.1,10},
PlotRange -> {{-1,5},{0.9,1}},AspectRatio->1]

Show[
allcplot
,Plot[ Table[c[2/p],{p,3,11,2}],{h,-1,5}]
, Plot[{1,25,1 + 24 h},{h,-15,15},PlotStyle->Dashed]
,PlotRange -> {{-1,5},{-20,2}}]

Show[
allcplot
,Plot[ Table[c[1/p],{p,2,11,1}],{h,-1,5}],PlotRange -> {{-1,5},{-20,2}}]

lines = Table[{(1/2)-((r^2)/2)-((Cos[a])/2)+(1/2) (r^2) Cos[a],13-12 Cos[a]},{r,1,50}];

lineplot = ParametricPlot[ %,{a,0,Pi}];

Show[allcplot,lineplot]

Constraints in the Lee-Yang model from null-vectors

h = -1/5; cc = -22/5 ;

initialise

sethwt

This calculates  <(Subscript[ϕ, 3 ])| (Subscript[ϕ, 2])(1)  (Subscript[ϕ, 3])(z) |ψ >
We can look at the three point function if we set (Subscript[h, 2 ])-> 0 and z-> 1

Simplify[diffeq[2]  /. h2 -> 0 ]
Factor[Simplify[%  /. y -> Function[z,z^(h3-h1+1/5)]] /. z -> 1]
equationslevel2 = ( %==0)

Simplify[diffeq[3]  /. h2 -> 0  ]
Factor[Simplify[% /. y -> Function[z,z^(h3-h1+1/5)]] /. z -> 1]
equationslevel3 = ( %==0)

ContourPlot[ {(1/25) (-2+25 (h1^2)-5 h3+25 (h3^2)-5 h1 (1+10 h3)),(1/5) (5 (h1^3)+h3+4 (h3^2)-5 (h3^3)-(h1^2) (4+15 h3)+h1 (-1+15 (h3^2)))},{h1,-.4,0.2},{h3,-.4,0.2}]

Solve[ {equationslevel2,equationslevel3}]