MMM_04.txt

In[1]:= (* Needs["VectorAnalysis`"]; *)
Needs["VariationalMethods`"];

generateMatrix[name_, n_, m_]:=Table[ToString[name] <> "[" <>ToString[ii - 1] <> ", " <>ToString[jj - 1] <> "]",{ii, 1, n}, {jj, 1, m}];
generateMatrix[name_, n_]:=generateMatrix[name, n, n];
generateMatrix[name_]:=generateMatrix[name, 3];

generateVector[name_, n_, m_]:=Table[ToString[name] <>ToString[ii - 1] <> ".["  <>ToString[jj - 1] <> "]",{ii, 1, n}, {jj, 1, m}];

complexFuncName="createComplex";
matrix3x3Name="ComplexMatrix3x3";
matrix4x4Name="ComplexMatrix4x4";
createName="create";
createFromReName="fromRe";
createFromImName="fromIm";

toFloat[s_?StringQ]:=If[StringCount[s,"."]==0,s<>".0",s];

toFSharpMatrix[m_?MatrixQ]:=
Module[{s,nn,mm,ii,jj,f,tolerance},
tolerance=10^-8;
s="\n[\n";
nn=Length[m];
mm=Length[m[[1]]];

Do[
(
If[mm>1,s=s<>"    [ "];

Do[
(
f=StringReplace[ToString[InputForm[m[[ii,jj]]]],"\n"->""];
s=s<>f<>If[jj<mm,"; ",""];
),{jj,1,mm}
];

If[mm>1,s=s<>" ]\n"];
),{ii,1,nn}
];

s=StringReplace[s<>"]\n","\""->""];
s=s<>"\n";
Return[s];
];

M[epsilon:{{_,_,_},{_,_,_},{_,_,_}},mu:{{_,_,_},{_,_,_},{_,_,_}},ro:{{_,_,_},{_,_,_},{_,_,_}},rotr:{{_,_,_},{_,_,_},{_,_,_}}]:=BlockMatrix[{{epsilon,ro},{rotr,mu}}];

BlockMatrix[{{AAA_,BBB_},{CCC_,DDD_}}]:=ArrayFlatten[{{AAA,BBB},{CCC,DDD}}];

OFULL[f_][x_,y_,z_]:=Transpose[{Join[Curl[{(f[x,y,z])[[4,1]],(f[x,y,z])[[5,1]],(f[x,y,z])[[6,1]]},{x,y,z}],-Curl[{(f[x,y,z])[[1,1]],(f[x,y,z])[[2,1]],(f[x,y,z])[[3,1]]},{x,y,z}]]}];


(* ==============================================*)
MMM[eps_,mu_,ro_,rotr_,lambda_,fita_,n1_]:=Module[{kx,MaxwellEq,EH,lst3,lst6,sol,xxxNoz,xNoz4,EH4,EH4f,EH4d,EH4df,Dl,Dld,MM,MMd,MMdInv,retval},

Print["MMM"];
Print["eps = ",MatrixForm[N[eps]]];
Print["mu = ",MatrixForm[N[mu]]];
Print["ro = ",MatrixForm[ro]];
Print["rotr = ",MatrixForm[rotr]];

EH[x_,y_,z_]:=Exp[I*kx*x]*{{Ex0[z]},{Ey0[z]},{Ez0[z]},{Hx0[z]},{Hy0[z]},{Hz0[z]}};

Print["EH[x,y,z] = ",EH[x,y,z]];

MaxwellEq=(((I*(Exp[(-I*kx*x)]*((2*Pi*I/lambda)*M[eps,mu,ro,rotr].EH[x,y,z]+OFULL[EH][x,y,z])))))/.x->0;
(* kx=(2*Pi/lambda)*n1*Sin[fita]; *)
kx=(2*Pi/lambda)*n1SinFita;

Print["kx = ",N[kx]];Print["n1 = ",N[n1]];
Print["fita = ",N[fita]];
Print["MaxwellEq = ",MatrixForm[Simplify[MaxwellEq]]];

lst3=CoefficientList[MaxwellEq[[3,1]],{Ez0[z],Hz0[z]}];
lst6=CoefficientList[MaxwellEq[[6,1]],{Ez0[z],Hz0[z]}];
(*
If[UseQuietSolveValue,(sol=(Quiet[Solve[{MaxwellEq[[3,1]]0,MaxwellEq[[6,1]]0},{Ez0[z],Hz0[z]}]])[[1]];),(sol=(Solve[{MaxwellEq[[3,1]]0,MaxwellEq[[6,1]]0},{Ez0[z],Hz0[z]}])[[1]];)];
*)
sol=(Solve[{MaxwellEq[[3,1]]==0,MaxwellEq[[6,1]]==0},{Ez0[z],Hz0[z]}])[[1]];
xxxNoz=MaxwellEq/.sol;
xNoz4={{xxxNoz[[1,1]]},{xxxNoz[[2,1]]},{xxxNoz[[4,1]]},{xxxNoz[[5,1]]}};
EH4={{Ex0[z]},{Hy0[z]},{Ey0[z]},{Hx0[z]}};
EH4f={{Ex0[z]},{Hy0[z]},{Ey0[z]},{-Hx0[z]}};
EH4d={{Derivative[1][Ex0][z]},{Derivative[1][Hy0][z]},{Derivative[1][Ey0][z]},{Derivative[1][Hx0][z]}};
EH4df={{Derivative[1][Ex0][z]},{Derivative[1][Hy0][z]},{Derivative[1][Ey0][z]},{-Derivative[1][Hx0][z]}};
Dl[q_,s_]:=(-Coefficient[xNoz4[[q,1]],EH4[[s,1]]]*If[s==4,-1,1]);
Dld[q_,s_]:=(Coefficient[xNoz4[[q,1]],EH4d[[s,1]]]*If[s==4,-1,1]);
MM=Simplify[Table[Dl[i,j],{i,4},{j,4}]];
MMd=Table[Dld[i,j],{i,4},{j,4}];
Print["MM = ", MM // MatrixForm] ;
Print["MMd = ", MMd // MatrixForm] ;

MMdInv=Inverse[MMd];
Print["MMdInv = ", MMdInv // MatrixForm] ;

retval=FullSimplify[{(MMdInv.MM)/(2*Pi*I/lambda),sol}];
(* retval=FullSimplify[(MMdInv . MM)/(2*Pi*I/lambda)]; *)
Print["MMM retval[[1]] = ",retval[[1]] // InputForm ];
Print["MMM retval[[2]] = ",retval[[2]] // InputForm ];
Return[retval];
];

(* ==============================================*)
SolutionNewBase[]:=
Module[{eps,ro,rotr,mu,eps2,ro2,rotr2,mu2,lambda,fita,n1,n2,beta,gamm,kx,Film,ehIncd,opts,calcBS,len,EHIv,EHRv,EHTv,ehr1,ehr2,eht1,eht2,ssss,delta,calcDlt,FilmLayer,EI,EP,ES,MMM1,MMM2,h2,PPPm,PPPmm,outPPPm,useSolve,sss1,sss2,cf,b,cmf,sol,varLst,free,coeffTbl,freeTbl,s1,s2,s3,s4,EGVal1Hlp,EGVal2Hlp,ODEG1hlp,EGVal1HlpIm,EGVal1HlpRe,ODEG2hlp,EGVal2HlpIm,EGVal2HlpRe,EGUnit,egs1,egs2,egvf1,egvf2,EHTFullEG1,EHTFullEG2,EHTFullEG3,EHTFullEG4,tmp,egvf1Tr,egvf2Tr,useNumEV,pdi,pdil,ehirule,ehrrule,ehtrule,EGVal1,EGVec1,ODEG1,EGVal1Up,EGVal1Dn,EGVec1Up,EGVec1Dn,EGVal2,EGVec2,ODEG2,EGVal2Up,EGVec2Up,EGVal2Dn,EGVec2Dn,ehi1,ehi2,ehisol,cfm,EHIcoeff,EHI,PPPv,coeff,ehrtsol,EHRcoeff,EHR,EHTcoeff,EHT,ehrule,EHIFull,EHTFull,EHRFull,retval, s },

kx=(2*Pi/lambda)*n1*Sin[fita];

MMM1=generateMatrix["m1.", 4];
MMM2=generateMatrix["m2.", 4];
(* 
Print["SolutionNewBase::MMM1 = ",MMM1 // MatrixForm];
Print["SolutionNewBase::MMM2 = ",MMM2 // MatrixForm];
*)

ehirule={Ez0[z] -> e1x * Ex0[z] + e1y * Ey0[z]+h1x * Hx0[z] + h1y * Hy0[z]};
ehrrule=ehirule;
ehtrule={Ez0[z] -> e2x * Ex0[z] + e2y * Ey0[z]+h2x * Hx0[z] + h2y * Hy0[z]};

(*
If[useNumEV===True,EGVal1=Eigenvalues[N[MMM1[[1]]]],EGVal1=Eigenvalues[MMM1[[1]]]];
If[pdiTrue&&pdil≥PCDILEVELDETAILED,Print["SolutionNewBase::EGVal1 = ",MatrixForm[EGVal1]];];
If[useNumEV===True,EGVec1=Eigenvectors[N[MMM1[[1]]]],EGVec1=Eigenvectors[MMM1[[1]]]];
If[pdiTrue&&pdil≥PCDILEVELALL,Print["SolutionNewBase::MMM1[[1]] = ",MatrixForm[Chop[N[MMM1[[1]]]]]];
Print["SolutionNewBase::EGVec1 = ",MatrixForm[EGVec1]];];

EGVal1Hlp=Chop[N[EGVal1]];
ODEG1=EGGetOrder[EGVal1Hlp,EGVec1,ehirule,opts];
EGVal1Up=EGVal1[[ODEG1[[1]]]];
EGVal1Dn:=EGVal1[[ODEG1[[2]]]];
EGVec1Up:=Transpose[(EGVec1[[ODEG1[[1]]]])];
EGVec1Dn:=Transpose[(EGVec1[[ODEG1[[2]]]])];
If[pdiTrue&&pdil≥PCDILEVELALL,Print["SolutionNewBase::ODEG1 = ",ODEG1];
Print["SolutionNewBase::ODEG1 Ended..."];];
If[useNumEV===True,EGVal2=Eigenvalues[N[MMM2[[1]]]],EGVal2=Eigenvalues[MMM2[[1]]]];

If[useNumEV===True,EGVec2=Eigenvectors[N[MMM2[[1]]],QuarticsTrue],EGVec2=Eigenvectors[MMM2[[1]],QuarticsTrue]];
EGVal2Hlp=Chop[N[EGVal2]];
ODEG2=EGGetOrder[EGVal2Hlp,EGVec2,ehtrule,opts];
EGVal2Up=EGVal2[[ODEG2[[1]]]];
EGVal2Dn=EGVal2[[ODEG2[[2]]]];
EGVec2Up=Transpose[(EGVec2[[ODEG2[[1]]]])];
EGVec2Dn=Transpose[(EGVec2[[ODEG2[[2]]]])];
*)

EGVec1Dn=Transpose[generateVector["b1.down.e", 2, 4]];
EGVec1Up=Transpose[generateVector["b1.up.e", 2, 4]];

EGVec2Up=Transpose[generateVector["b2.down.e", 2, 4]];
EGVec2Dn=Transpose[generateVector["b2.up.e", 2, 4]];

(*
Print["EGVec1Dn = ", EGVec1Dn // MatrixForm];
Print["EGVec1Up = ", EGVec1Up // MatrixForm];
Print["EGVec2Up = ", EGVec2Up // MatrixForm];
Print["EGVec2Dn = ", EGVec2Dn // MatrixForm];
*)


(* EGVec1Dn=Transpose[{{e111, e112, e113, e114}}]; *)
(* EGVec1Dn=Transpose[{{e111, e112},{e121, e122}}]; *)

(* ehIncd={eXi, eYi, eZi, hXi, hYi, hZi}; *)
ehIncd={"i.e.x", "i.e.y", "i.e.z", "i.h.x","i.h.y", "i.h.z"};

EI={{ehIncd[[1]]},{ehIncd[[5]]},{ehIncd[[2]]},{-ehIncd[[4]]}};
(* Print["EI = ", EI // MatrixForm]; *)

Clear[ehi1,ehi2,ssss];
ssss=(EGVec1Dn . {{ehi1},{ehi2}}-EI);
(* Print["ssss = ", ssss // MatrixForm]; *)

sss1=ssss[[1,1]];
sss2=ssss[[3,1]];
cf={{Coefficient[sss1,ehi1],Coefficient[sss1,ehi2]},{Coefficient[sss2,ehi1],Coefficient[sss2,ehi2]}};
b=-{{sss1/.ehi1->0/.ehi2->0},{sss2/.ehi1->0/.ehi2->0}};
(*
Print["sss1 = ", sss1];
Print["sss2 = ", sss2];
Print["cf = ", cf // MatrixForm];
Print["b = ", b // MatrixForm];
*)
cfm=Inverse[cf];
(* Print["cfm = ", cfm // MatrixForm]; *)

sol=Simplify[cfm.b];
(* Print["sol = ", sol // MatrixForm]; *)

ehisol={ehi1->sol[[1,1]],ehi2->sol[[2,1]]};
(* Print["ehisol = ", ehisol // MatrixForm]; *)

EHIcoeff=(Clear[ehi1,ehi2];({{ehi1},{ehi2}}/.ehisol));
EHI=EGVec1Dn.EHIcoeff;
(*
Print["EHIcoeff = ", EHIcoeff // MatrixForm];
Print["EHI = ", EHI // MatrixForm];
*)
(*Print["PPPv calc started..."];*)

(* PPPv=If[calcBS===True,II,PPPFull[Film,lambda,fita,n1],II]; *)
PPPv=generateMatrix["p.", 4];

EHIv=EHI;
EHRv=EGVec1Up.{{ehr1},{ehr2}};
EHTv=EGVec2Dn.{{eht1},{eht2}};
ssss=(PPPv.(EHIv+EHRv)-EHTv);
(*
Print["EHIv = ", EHIv // MatrixForm];
Print["EHRv = ", EHRv // MatrixForm];
Print["EHTv = ", EHTv // MatrixForm];
Print["ssss = ", ssss // MatrixForm];
*)
varLst={ehr1,ehr2,eht1,eht2};
coeff[q_,s_]:=Coefficient[ssss[[q,1]],varLst[[s]]];
free[q_]:=(-ssss[[q,1]]/.ehr1->0/.ehr2->0/.eht1->0/.eht2->0);
coeffTbl=FullSimplify[Table[coeff[i,j],{i,4},{j,4}]];

(* Print["varLst = ", varLst // MatrixForm]; *)
Print["coeffTbl = ", coeffTbl // toFSharpMatrix];

freeTbl=FullSimplify[Table[free[i],{i,4},{j,1}]];
cfm=Inverse[coeffTbl];


(* Print["cfm = ", cfm // toFSharpMatrix]; *)
Print["freeTbl = ", freeTbl // toFSharpMatrix];


freeTbl=generateMatrix["ft.", 4,1];
cfm=generateMatrix["cfm.", 4];

Print["cfm (subst) = ", cfm // toFSharpMatrix];
Print["freeTbl (subst) = ", freeTbl // toFSharpMatrix];


s1=ssss[[1,1]];
s2=ssss[[2,1]];
s3=ssss[[3,1]];
s4=ssss[[4,1]];

sol=cfm.freeTbl;
Print["sol = ", sol // MatrixForm];

sol=generateMatrix["sol.", 4,1];
Print["sol = ", sol // MatrixForm];

ehrtsol={ehr1->sol[[1,1]],ehr2->sol[[2,1]],eht1->sol[[3,1]],eht2->sol[[4,1]]};
Print["ehrtsol = ", ehrtsol // MatrixForm];


EHRcoeff=(Clear[ehr1,ehr2];({{ehr1},{ehr2}}/.ehrtsol));
EHTcoeff=(Clear[eht1,eht2];({{eht1},{eht2}}/.ehrtsol));
EHR=EGVec1Up.EHRcoeff;
EHT=EGVec2Dn.EHTcoeff;
Print["EHR = ", EHR // toFSharpMatrix];
Print["EHT = ", EHT // toFSharpMatrix];

Abort[];

ehrule={Ex0[z]->EHI[[1,1]],Ey0[z]->EHI[[3,1]],Hx0[z]->-EHI[[4,1]],Hy0[z]->EHI[[2,1]]};
EHIFull={EHI[[1,1]],EHI[[3,1]],((Ez0[z]/.ehirule)/.ehrule),-EHI[[4,1]],EHI[[2,1]],((Hz0[z]/.ehirule)/.ehrule),True};
ehrule=({Ex0[z]->EHT[[1,1]],Ey0[z]->EHT[[3,1]],Hx0[z]->-EHT[[4,1]],Hy0[z]->EHT[[2,1]]});
EHTFull={EHT[[1,1]],EHT[[3,1]],((Ez0[z]/.ehtrule)/.ehrule),-EHT[[4,1]],EHT[[2,1]],((Hz0[z]/.ehtrule)/.ehrule),True};
ehrule=({Ex0[z]->EHR[[1,1]],Ey0[z]->EHR[[3,1]],Hx0[z]->-EHR[[4,1]],Hy0[z]->EHR[[2,1]]});
EHRFull={EHR[[1,1]],EHR[[3,1]],((Ez0[z]/.ehrrule)/.ehrule),-EHR[[4,1]],EHR[[2,1]],((Hz0[z]/.ehrrule)/.ehrule),False};
(*PCDILEVELALL;PCDILEVELDETAILED;PCDILEVELMEDIUM;PCDILEVELSHORT;*)If[pdi==True&&pdil>=PCDILEVELALL,Print["EHIFull = ",Chop[N[EHIFull]]];
Print["EHTFull = ",Chop[N[EHTFull]]];
Print["EHRFull = ",Chop[N[EHRFull]]];];
If[calcDlt===True&&len>0,FilmLayer=Film[[1]];
eps=FilmLayerEpsilon[FilmLayer];
mu=FilmLayerMu[FilmLayer];
ro=FilmLayerRo[FilmLayer];
rotr=FilmLayerRoT[FilmLayer];
delta=(MMM[eps,mu,ro,rotr,lambda,fita,n1])[[1]];];
egvf1=JoinRight[EGVec1Up,EGVec1Dn];
egvf2=JoinRight[EGVec2Up,EGVec2Dn];
egvf1Tr=Transpose[egvf1];
egvf2Tr=Transpose[egvf2];
ehrule=({Ex0[z]->EHR[[1,1]],Ey0[z]->EHR[[3,1]],Hx0[z]->-EHR[[4,1]],Hy0[z]->EHR[[2,1]]});
tmp={egvf2[[1,1]],egvf2[[3,1]],(Ez0[z]),-egvf2[[4,1]],egvf2[[2,1]],(Hz0[z]),True};
(*Print["tmp = ",tmp];*)EHTFullEG1=GetEHFull[egvf2Tr[[1]],ehtrule,True];
EHTFullEG2=GetEHFull[egvf2Tr[[2]],ehtrule,True];
EHTFullEG3=GetEHFull[egvf2Tr[[3]],ehtrule,True];
EHTFullEG4=GetEHFull[egvf2Tr[[4]],ehtrule,True];
(*EHTFullEG1={egvf2[[1,1]],egvf2[[3,1]],((Ez0[z]/.ehtrule)/.ehrule),-egvf2[[4,1]],egvf2[[2,1]],((Hz0[z]/.ehtrule)/.ehrule),True};
EHTFullEG2={egvf2[[1,2]],egvf2[[3,2]],((Ez0[z]/.ehtrule)/.ehrule),-egvf2[[4,2]],egvf2[[2,2]],((Hz0[z]/.ehtrule)/.ehrule),True};
EHTFullEG3={egvf2[[1,3]],egvf2[[3,3]],((Ez0[z]/.ehtrule)/.ehrule),-egvf2[[4,3]],egvf2[[2,3]],((Hz0[z]/.ehtrule)/.ehrule),True};
EHTFullEG4={egvf2[[1,4]],egvf2[[3,4]],((Ez0[z]/.ehtrule)/.ehrule),-egvf2[[4,4]],egvf2[[2,4]],((Hz0[z]/.ehtrule)/.ehrule),True};*)egs1={Flatten[{EGVal1Up,EGVal1Dn}],egvf1};
egs2={Flatten[{EGVal2Up,EGVal2Dn}],egvf2,EHTFullEG1,EHTFullEG2,EHTFullEG3,EHTFullEG4};
(*Print["fita = ",Chop[N[fita/Degree]]];
Print["egs1 = ",Chop[N[egs1]]];Print["egs2 = ",Chop[N[egs2]]];*)retval={EHIFull,EHRFull,EHTFull,PPPv,delta,Media,IncidentLight,opts,MMM1[[1]],MMM2[[1]],coeffTbl,freeTbl,egs1,egs2};
(*PCDILEVELALL;PCDILEVELDETAILED;PCDILEVELMEDIUM;PCDILEVELSHORT;*)If[pdi==True&&pdil>=PCDILEVELALL,Print["SolutionNewBase::retval = ",retval];];
If[pdi==True,Print["SolutionNewBase::end ================================================="];
Print["   "];];
Return[retval];];

(*
eps = {{e11,e12,e13},{e21, e22, e23},{e31, e32, e33}};
mu = {{m11,m12,m13},{m21, m22, m23},{m31, m32, m33}};
rho = {{r11,r12,r13},{r21, r22, r23},{r31, r23, e33}};
rhoT =  {{r11,r12,r13},{r21, r22, r23},{r31, r23, e33}};
*)
eps = generateMatrix["o.eps."];
mu =  generateMatrix["o.mu."];
rho =  generateMatrix["o.rho."];
rhoT =  generateMatrix["o.rhoT."];

eps // MatrixForm
(* retVal=MMM[eps, mu, rho, rhoT, lambda,fita,n1]; *)

Print["============================="];

SolutionNewBase[];

Out[23]//MatrixForm= (o.eps.[0, 0]	o.eps.[0, 1]	o.eps.[0, 2]
o.eps.[1, 0]	o.eps.[1, 1]	o.eps.[1, 2]
o.eps.[2, 0]	o.eps.[2, 1]	o.eps.[2, 2]

)
During evaluation of In[1]:= =============================
During evaluation of In[1]:= coeffTbl = 
[
    [ b1.up.e0.[0]*p.[0, 0] + b1.up.e0.[1]*p.[0, 1] + b1.up.e0.[2]*p.[0, 2] + b1.up.e0.[3]*p.[0, 3]; b1.up.e1.[0]*p.[0, 0] + b1.up.e1.[1]*p.[0, 1] + b1.up.e1.[2]*p.[0, 2] + b1.up.e1.[3]*p.[0, 3]; -b2.up.e0.[0]; -b2.up.e1.[0] ]
    [ b1.up.e0.[0]*p.[1, 0] + b1.up.e0.[1]*p.[1, 1] + b1.up.e0.[2]*p.[1, 2] + b1.up.e0.[3]*p.[1, 3]; b1.up.e1.[0]*p.[1, 0] + b1.up.e1.[1]*p.[1, 1] + b1.up.e1.[2]*p.[1, 2] + b1.up.e1.[3]*p.[1, 3]; -b2.up.e0.[1]; -b2.up.e1.[1] ]
    [ b1.up.e0.[0]*p.[2, 0] + b1.up.e0.[1]*p.[2, 1] + b1.up.e0.[2]*p.[2, 2] + b1.up.e0.[3]*p.[2, 3]; b1.up.e1.[0]*p.[2, 0] + b1.up.e1.[1]*p.[2, 1] + b1.up.e1.[2]*p.[2, 2] + b1.up.e1.[3]*p.[2, 3]; -b2.up.e0.[2]; -b2.up.e1.[2] ]
    [ b1.up.e0.[0]*p.[3, 0] + b1.up.e0.[1]*p.[3, 1] + b1.up.e0.[2]*p.[3, 2] + b1.up.e0.[3]*p.[3, 3]; b1.up.e1.[0]*p.[3, 0] + b1.up.e1.[1]*p.[3, 1] + b1.up.e1.[2]*p.[3, 2] + b1.up.e1.[3]*p.[3, 3]; -b2.up.e0.[3]; -b2.up.e1.[3] ]
]


During evaluation of In[1]:= freeTbl = 
[
((b1.down.e1.[2]*i.e.x - b1.down.e1.[0]*i.e.y)*(b1.down.e0.[1]*p.[0, 1] + b1.down.e0.[3]*p.[0, 3]) - b1.down.e0.[2]*(b1.down.e1.[0]*i.e.x*p.[0, 0] + b1.down.e1.[1]*i.e.x*p.[0, 1] + b1.down.e1.[0]*i.e.y*p.[0, 2] + b1.down.e1.[3]*i.e.x*p.[0, 3]) + b1.down.e0.[0]*(b1.down.e1.[2]*i.e.x*p.[0, 0] + b1.down.e1.[1]*i.e.y*p.[0, 1] + b1.down.e1.[2]*i.e.y*p.[0, 2] + b1.down.e1.[3]*i.e.y*p.[0, 3]))/(b1.down.e0.[2]*b1.down.e1.[0] - b1.down.e0.[0]*b1.down.e1.[2])((b1.down.e1.[2]*i.e.x - b1.down.e1.[0]*i.e.y)*(b1.down.e0.[1]*p.[1, 1] + b1.down.e0.[3]*p.[1, 3]) - b1.down.e0.[2]*(b1.down.e1.[0]*i.e.x*p.[1, 0] + b1.down.e1.[1]*i.e.x*p.[1, 1] + b1.down.e1.[0]*i.e.y*p.[1, 2] + b1.down.e1.[3]*i.e.x*p.[1, 3]) + b1.down.e0.[0]*(b1.down.e1.[2]*i.e.x*p.[1, 0] + b1.down.e1.[1]*i.e.y*p.[1, 1] + b1.down.e1.[2]*i.e.y*p.[1, 2] + b1.down.e1.[3]*i.e.y*p.[1, 3]))/(b1.down.e0.[2]*b1.down.e1.[0] - b1.down.e0.[0]*b1.down.e1.[2])((b1.down.e1.[2]*i.e.x - b1.down.e1.[0]*i.e.y)*(b1.down.e0.[1]*p.[2, 1] + b1.down.e0.[3]*p.[2, 3]) - b1.down.e0.[2]*(b1.down.e1.[0]*i.e.x*p.[2, 0] + b1.down.e1.[1]*i.e.x*p.[2, 1] + b1.down.e1.[0]*i.e.y*p.[2, 2] + b1.down.e1.[3]*i.e.x*p.[2, 3]) + b1.down.e0.[0]*(b1.down.e1.[2]*i.e.x*p.[2, 0] + b1.down.e1.[1]*i.e.y*p.[2, 1] + b1.down.e1.[2]*i.e.y*p.[2, 2] + b1.down.e1.[3]*i.e.y*p.[2, 3]))/(b1.down.e0.[2]*b1.down.e1.[0] - b1.down.e0.[0]*b1.down.e1.[2])((b1.down.e1.[2]*i.e.x - b1.down.e1.[0]*i.e.y)*(b1.down.e0.[1]*p.[3, 1] + b1.down.e0.[3]*p.[3, 3]) - b1.down.e0.[2]*(b1.down.e1.[0]*i.e.x*p.[3, 0] + b1.down.e1.[1]*i.e.x*p.[3, 1] + b1.down.e1.[0]*i.e.y*p.[3, 2] + b1.down.e1.[3]*i.e.x*p.[3, 3]) + b1.down.e0.[0]*(b1.down.e1.[2]*i.e.x*p.[3, 0] + b1.down.e1.[1]*i.e.y*p.[3, 1] + b1.down.e1.[2]*i.e.y*p.[3, 2] + b1.down.e1.[3]*i.e.y*p.[3, 3]))/(b1.down.e0.[2]*b1.down.e1.[0] - b1.down.e0.[0]*b1.down.e1.[2])]


During evaluation of In[1]:= cfm (subst) = 
[
    [ cfm.[0, 0]; cfm.[0, 1]; cfm.[0, 2]; cfm.[0, 3] ]
    [ cfm.[1, 0]; cfm.[1, 1]; cfm.[1, 2]; cfm.[1, 3] ]
    [ cfm.[2, 0]; cfm.[2, 1]; cfm.[2, 2]; cfm.[2, 3] ]
    [ cfm.[3, 0]; cfm.[3, 1]; cfm.[3, 2]; cfm.[3, 3] ]
]


During evaluation of In[1]:= freeTbl (subst) = 
[
ft.[0, 0]ft.[1, 0]ft.[2, 0]ft.[3, 0]]


During evaluation of In[1]:= sol = (cfm.[0, 0] ft.[0, 0]+cfm.[0, 1] ft.[1, 0]+cfm.[0, 2] ft.[2, 0]+cfm.[0, 3] ft.[3, 0]
cfm.[1, 0] ft.[0, 0]+cfm.[1, 1] ft.[1, 0]+cfm.[1, 2] ft.[2, 0]+cfm.[1, 3] ft.[3, 0]
cfm.[2, 0] ft.[0, 0]+cfm.[2, 1] ft.[1, 0]+cfm.[2, 2] ft.[2, 0]+cfm.[2, 3] ft.[3, 0]
cfm.[3, 0] ft.[0, 0]+cfm.[3, 1] ft.[1, 0]+cfm.[3, 2] ft.[2, 0]+cfm.[3, 3] ft.[3, 0]

)
During evaluation of In[1]:= sol = (sol.[0, 0]
sol.[1, 0]
sol.[2, 0]
sol.[3, 0]

)
During evaluation of In[1]:= ehrtsol = (ehr1$2164->sol.[0, 0]
ehr2$2164->sol.[1, 0]
eht1$2164->sol.[2, 0]
eht2$2164->sol.[3, 0]

)
During evaluation of In[1]:= EHR = 
[
b1.up.e0.[0]*sol.[0, 0] + b1.up.e1.[0]*sol.[1, 0]b1.up.e0.[1]*sol.[0, 0] + b1.up.e1.[1]*sol.[1, 0]b1.up.e0.[2]*sol.[0, 0] + b1.up.e1.[2]*sol.[1, 0]b1.up.e0.[3]*sol.[0, 0] + b1.up.e1.[3]*sol.[1, 0]]


During evaluation of In[1]:= EHT = 
[
b2.up.e0.[0]*sol.[2, 0] + b2.up.e1.[0]*sol.[3, 0]b2.up.e0.[1]*sol.[2, 0] + b2.up.e1.[1]*sol.[3, 0]b2.up.e0.[2]*sol.[2, 0] + b2.up.e1.[2]*sol.[3, 0]b2.up.e0.[3]*sol.[2, 0] + b2.up.e1.[3]*sol.[3, 0]]


Out[25]= $Aborted