#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include"comple.h"

 void main(COMPLEX A[4][4],double d,double ramda,COMPLEX ep1,COMPLEX ep2)
{
 COMPLEX epp, enn, dp, dn, sp, sn, cp, cn;
 COMPLEX spd, snd, spm, snm, acpcn, scpcn, aspdsnd, sspdsnd;
 COMPLEX aspmsnm, sspmsnm;

 epp=csqrt(cadd(ep1, ep2));
 enn=csqrt(csub(ep1, ep2));

 dp=cmul_r(epp,(2.0*PAI*d/ramda));
 dn=cmul_r(enn,(2.0*PAI*d/ramda));

 sp=zsin(dp);
 sn=zsin(dn);
 cp=zcos(dp);
 cn=zcos(dn);

 spd=cdiv(sp, epp);
 snd=cdiv(sn, enn);
 spm=cmul(sp, epp);
 snm=cmul(sn, enn);

 acpcn=cadd(cp, cn);
 scpcn=csub(cp, cn);
 aspdsnd=cadd(spd, snd);
 sspdsnd=csub(spd, snd);
 aspmsnm=cadd(spm, snm);
 sspmsnm=csub(spm, snm);

 clear_mat(&A[0][0],4,4);

 A[0][0]=cmul_r(acpcn,0.5);
 A[0][1]=cmul(scpcn,cmplx(0.0,0.5));
 A[0][2]=cmul_r(sspdsnd,0.5);
 A[0][3]=cmul(aspdsnd,cmplx(0.0,0.5));

 A[1][0]=cmul(scpcn,cmplx(0.0,-0.5));
 A[1][2]=cmul_r(acpcn,0.5);
 A[1][2]=cmul(aspdsnd, cmplx(0.0,-0.5));
 A[1][3]=cmul_r(sspdsnd,0.5);

 A[2][0]=cmul_r(sspmsnm,-0.5);
 A[2][1]=cmul(aspmsnm,cmplx(0.0,-0.5));
 A[2][2]=cmul_r(acpcn,0.5);
 A[2][3]=cmul(scpcn,cmplx(0.0,0.5));

 A[3][0]=cmul(aspmsnm,cmplx(0.0,0.5));
 A[3][1]=cmul_r(sspmsnm,-0.5);
 A[3][2]=cmul(scpcn,cmplx(0.0,-0.5));
 A[3][3]=cmul_r(acpcn,0.5);
 printf("yaho!");
 }


 

 #include<math.h>
#define CZERO cmplx(0.0,0.0)
#define CONE cmplx(1.0,0.0)
#define ZI cmplx(0.0,1.0)
#define PAI 3.1415926535e+00
#define cabs(z) sqrt(ab2((z)))


 COMPLEX cadd(COMPLEX a,COMPLEX b)
{
 COMPLEX z;
 z.re=a.re+b.re;
 z.im=a.im+b.im;
 return(z);
 }

COMPLEX csub(COMPLEX a,COMPLEX b)
{
 COMPLEX z;
 z.re=a.re-b.re;
 z.im=a.im-b.im;
 return(z);
 }

COMPLEX cmul(COMPLEX a,COMPLEX b)
{
 COMPLEX z;
 z.re=a.re*b.re-a.im*b.im;
 z.im=a.im*b.re+a.re*b.im;
 return(z);
 }
COMPLEX cdiv(COMPLEX a,COMPLEX b)
{
 double ab;
 COMPLEX z;
 ab=ab2(b);
 z=cmul(a,conj(b));
 z.re /= ab;
 z.im /=ab;
 return(z);
 }

COMPLEX cmul_r(COMPLEX z,double r)
{
 COMPLEX zr;
 zr.re=z.re*r;
 zr.im=z.im*r;
 return(zr);
 }
COMPLEX csqrt(COMPLEX z)
{
 COMPLEX zr;
 double fai,r;
 if(z.re==0.0){
  if(z.im==0.0){
   fai=0.0;
   }
  else if(z.im>0.0){
   fai=PAI/4.0;
   }
  else{
   fai=PAI*0.75;
   }
  }
 else {
  fai=(atan(z.im/z.re))/2.0;
  if(z.re<0.0){
   fai+=PAI/2.0;
   }
  if((z.re>0.0)&&(z.im<0.0)){
   fai+=PAI;
   }
  }
 zr.re=(r=sqrt(cabs(z)))*cos(fai);
 zr.im=r*sin(fai);
 return(zr);
 }

COMPLEX zsin(COMPLEX za)
{
 COMPLEX ep,en;
 ep=zexp(cmul(cmplx(0.0,1.0),za));
 en=zexp(cmul(cmplx(0.0,-1.0),za));
 return(cdiv(csub(ep,en),cmplx(0.0,2.0)));
 }

COMPLEX zcos(COMPLEX za)
{
 COMPLEX ep,en;
 ep=zexp(cmul(cmplx(0.0,1.0),za));
 en=zexp(cmul(cmplx(0.0,-1.0),za));
 return(cdiv(cadd(ep,en),cmplx(2.0,0.0)));
 }

 COMPLEX cmplx(double a,double b)
{
 COMPLEX z;
 z.re=a;
 z.im=b;
 return(z);
 }

 double ab2(COMPLEX z)
{
 return(z.re*z.re+z.im*z.im);
        }
COMPLEX conj(COMPLEX z)
{
 COMPLEX a;
 a.re=z.re;
 a.im=-z.im;
 return(a);
 }
void clear_mat(COMPLEX *a,int m,int n)
{
 COMPLEX *paend;
 for(paend=a+m*n;a<paend;)
  *a++=CZERO;
 }
COMPLEX zexp(COMPLEX z)
{
	COMPLEX zz;
	zz=cmul(cmplx(exp(z.re),0.0),cmplx(cos(z.im),sin(z.im)));
 return(zz);
 }