omega=2*pi/0.1492 R0=7.112 index=0.108 b=1 eta=R0/(1-index) Ap=0.045 a1=1; a2=1;a3=1. dw=0.0002*omega f(x,n,off)=cos((omega+dw*n)*(x+off)) gp(x,off)= a1*f(x,0,off)+a2*f(x,1,off)+a2*f(x,2,off)+a3*f(x,3,off)+f(x,4,off)+f(x,5,off)+f(x,6,off)+f(x,7,off)+f(x,8,off)+f(x,9,off)+f(x,10,off) gm(x,off)= a1*f(x,0,off)+a2*f(x,-1,off)+a2*f(x,-2,off)+a3*f(x,-3,off)+f(x,-4,off)+f(x,-5,off)+f(x,-6,off)+f(x,-7,off)+f(x,-8,off)+f(x,-9,off)+f(x,-10,off) h(x,a,sigma) = a*exp(-x**2/2/sigma**2) plot gp(x,0),gm(x,0) pause -1 g(x,off)=gp(x,off)+gm(x,off) plot g(x,0)+g(x,0.05)+g(x,-0.05) pause -1 f1(x,dmin,dmax)=(-eta/omega/b/x*(dmax+dmin)*(cos(omega*(1+b*dmax)*x)-cos(omega*(1+b*dmin)*x)))/(dmax-dmin) f2(x,d)=-(eta*sin(omega*x)*(1/(omega*b*x)**2*cos(omega*b*d*x)+d/(omega*b*x)*sin(omega*b*d*x)))/(dmax-dmin) f2both(x,dmin,dmax)=f2(x,dmax)-f2(x,dmin) f3(x,d)=(eta*cos(omega*x)*(1/(omega*b*x)**2*sin(omega*b*d*x)-d/(omega*b*x)*cos(omega*b*x*d)))/(dmax-dmin) f3both(x,dmin,dmax)=f3(x,dmax)-f3(x,dmin) dmax=Ap/2/eta dmin=-Ap/2/eta plot f1(x,dmin,dmax)+f2both(x,dmin,dmax)+f3both(x,dmin,dmax) ftot(x,dmin,dmax)=f1(x,dmin,dmax)+f2both(x,dmin,dmax)+f3both(x,dmin,dmax) pause -1 off=0.05 plot ftot(x+0.1,0,Ap/eta), amp*exp(-x**2/2/sigma**2)