#This file converts a given real Prandtl stress function phi to complex variable form.

#phi:=-B*nu*x^3/(6*(1+nu))+B*nu*a^2*x/(6*(1+nu))-B*a^2*x/2;

delphi:=simplify(diff(phi,x,x)+diff(phi,y,y));

#Note we use the alternative notation AA,BB,CC here to avoid confusion with A[0] etc.

AA:=coeff(delphi,y)*(1+nu)/nu;
BB:=-coeff(delphi,x)*(1+nu)/nu;
CC:=simplify(delphi-nu*(AA*y-BB*x)/(1+nu));

phi0:=-BB*(3*(1+2*nu)*y^2*x-(1-2*nu)*x^3)/(24*(1+nu))+AA*(3*(1+2*nu)*x^2*y-(1-2*nu)*y^3)/(24*(1+nu))+CC*(x^2+y^2)/4;

delphi0:=simplify(diff(phi0,x,x)+diff(phi0,y,y));

f:=simplify(phi-phi0);

read rtoch:

h:=2*I*g;
hb:=subs({z1=z2,z2=z1,I=-I},h);

read GH:

phi:=simplify(subs({A[0]=(AA-I*BB)/2,Ab[0]=(AA+I*BB)/2,B[0]=0},G[0,1])+z*(h+hb)-z*(z1*diff(h,z1)+z2*diff(hb,z2))/(4*(1-nu)));
psi:=simplify(subs({A[0]=(AA-I*BB)/2,Ab[0]=(AA+I*BB)/2,B[0]=0},H[0,1])+z*diff(hb,z2)/(2*(1-nu))+I*CC*chi[1,1]/(4*(1-nu)));
psib:=subs({z1=z2,z2=z1,I=-I},psi);

#The following lines, if uncommented, would calculate the stress components in both complex and real form.

#read compPN:

#szz:=simplify(szz);
#Phi:=simplify(Phi);
#Psi1:=simplify(Psi1);
#Theta:=simplify(Theta);

#Psi1:=expand(subs({z1=x+I*y,z2=x-I*y},Psi1));
#szx:=coeff(Psi1,I,0);
#szy:=coeff(Psi1,I,1);