function Y = Residual_linear(f_T1,f_T2,Z_select,Z)
global eps;global delv;
[Nx, Nv, n] = size(f_T1);N = Nx*Nv;
[Nz,~] = size(Z);
[sqrt_M_select,A_select] = fcn_M_A(Z_select);
[sqrt_M,A] = fcn_M_A(Z);


vdx_fn = repmat(vdx_fn_perd_1st(f_T1),1,1,1,Nz);
sig_11 = repmat(reshape(sig_1(Z),1,1,1,Nz),1,1,n,1);
sig_12 = repmat(reshape(sig_1(Z_select),1,1,n),1,1,1,Nz);
part_1 = (sig_11 - sig_12).*vdx_fn;

Pf_1 = Trig_Multi(A_select,f_T2./sqrt_M_select).*sqrt_M_select/(delv^2);
Pf_1 =  MultiD_times(Pf_1,sig_3(Z_select)/eps);
Pf_1 = repmat(Pf_1,1,1,1,Nz);
A = repmat(reshape(A,Nx,Nv,1,Nz),1,1,n,1);
sqrt_M = repmat(reshape(sqrt_M,Nx,Nv,1,Nz),1,1,n,1);
Pf_2 = Trig_Multi(A,f_T2./sqrt_M).*sqrt_M/(delv^2);
sig_32 = repmat(reshape(sig_3(Z),1,1,1,Nz),Nx,Nv,n,1);
Pf_2 =  Pf_2.*sig_32/eps;
part_2 = Pf_2 - Pf_1;

Y = reshape(part_1 - part_2,N,n,Nz); 

