function [taoEM]=TaoEM(y,sita)
  k=y^2;
  k1=sqrt(1+2*k*cos(sita)+k^2);
  sita1=pi+asin(k*sin(sita)/k1);
  s3=k*k1*cos(sita)*cos(sita1)-2*k*k1*cos(sita1-sita);
  s4=sqrt(k*k1*cos(sita1-sita))-(0.011-0.012*j)/2;
  taoEM=(1/2)*(s3/s4);