MATLAB建模仿真案例：26 模拟电磁场分布.zip

clear all BASTA=0; %serve a fermare c=299792458; hbar=6.5822e-16; pi=acos(-1); lambda=632.8e-9; omega=2*pi/lambda*c; fukso=complex(0,1); %setup del dialogo per i parametri: figure(1) clf axes('position',[.05 .3 .9 .68]) %PRISMA n_prisma=uicontrol('style','slider','min',0,'max',2,'value',1.505); set(n_prisma,'units','normalized','position',[.10 .16 .25 .025]); uicontrol('style','text','string','PRISM','units','normalized','position',[.0 .16 .07 .02]); %ORO n_oro=uicontrol('style','slider','min',0,'max',1,'value',.1726); set(n_oro,'units','normalized','position',[.10 .12 .25 .025]); k_oro=uicontrol('style','slider','min',0,'max',4,'value',3.4218); set(k_oro,'units','normalized','position',[.40 .12 .25 .025]); d_oro=uicontrol('style','slider','min',0,'max',100,'value',88); set(d_oro,'units','normalized','position',[.70 .12 .25 .025]); uicontrol('style','text','string','GOLD','units','normalized','position',[.0 .12 .07 .02]); %LAYER n_lyr=uicontrol('style','slider','min',0,'max',2,'value',1.7); set(n_lyr,'units','normalized','position',[.10 .08 .25 .025]); k_lyr=uicontrol('style','slider','min',0,'max',1,'value',0); set(k_lyr,'units','normalized','position',[.40 .08 .25 .025]); d_lyr=uicontrol('style','slider','min',0,'max',400,'value',20); set(d_lyr,'units','normalized','position',[.70 .08 .25 .025]); uicontrol('style','text','string','LAYER','units','normalized','position',[.0 .08 .07 .02]); % SUBPHASE CHOICE subfase=uicontrol('style','popup','string','air|water'); set(subfase,'units','normalized','position',[.13 .01 .1 .05]) uicontrol('style','text','string','SUBPHASE','units','normalized','position',[.0 .03 .1 .02]); %FILE CHOICE FL=dir('*.dat'); MS=' '; for j=1:length(FL); MS=[MS,'|',FL(j).name]; end fname=uicontrol('style','popup','string',MS); set(fname,'units','normalized','position',[.5 .01 .2 .05]) uicontrol('style','text','string','EXPERIMENTAL DATA FILE','units','normalized','position',[.3 .03 .2 .02]); %CAPTIONS uicontrol('style','text','string','n','units','normalized','position',[.2 .22 .02 .02]); uicontrol('style','text','string','k','units','normalized','position',[.5 .22 .02 .02]); uicontrol('style','text','string','d','units','normalized','position',[.8 .22 .02 .02]); %TERMINA CALCOLI basta=uicontrol('style','pushbutton','string','STOP','callback','BASTA=1;'); set(basta,'units','normalized','position',[.9 .01 .1 .05]); while BASTA==0 clear expdata REF TRA %se passo ad un dataset pi?corto ... %% file di dati speriementali da fittare fnumber=get(fname,'value'); if fnumber==1 expdata(:,1)=(30:.1:70)'; expdata(:,2)=1; FNAME='SIMULATION ONLY'; else FNAME=FL(fnumber-1).name; expdata=dlmread(FNAME,' ',15,1); end %strato 1: vetro del prisma en(1)=get(n_prisma,'value'); %parte reale indice rifrazione prisma (1.51 se vetro normale, 1.723 se SF10) ek(1)=0; %strato 2: layer di oro en(2)=get(n_oro,'value'); %parte reale indice rifr oro ek(2)=get(k_oro,'value'); %parte imag indice rifr oro d(2)=1e-9*get(d_oro,'value'); %spessore (nm) %strato 3: layer appiccicato en(3)=get(n_lyr,'value'); %parte reale indice rifrazione layer ek(3)=get(k_lyr,'value'); %parte imag indice rifrazione layer d(3)=1e-9*get(d_lyr,'value'); %spessore (nm) %strato 4: subfase sottostante if get(subfase,'value')==1 %AIR en(4)=1.00; ek(4)=0; else %WATER en(4)=1.33; ek(4)=0; end %range angolare su cui fare il conto THETA_ext_deg=expdata(:,1); %calcolo costanti dielettriche er=en(1)^2-ek(1)^2; ei=2*en(1)*ek(1); e(1)=complex(er,ei); er=en(2)^2-ek(2)^2; ei=2*en(2)*ek(2); e(2)=complex(er,ei); er=en(3)^2-ek(3)^2; ei=2*en(3)*ek(3); e(3)=complex(er,ei); er=en(4)^2-ek(4)^2; ei=2*en(4)*ek(4); e(4)=complex(er,ei); % --------- CALCOLI VERI ---------- THETA_ext=THETA_ext_deg/180*pi; THETA=pi/4+asin(1/en(1)*sin(THETA_ext-pi/4)); for jtheta=1:length(THETA); theta=THETA(jtheta); q1=sqrt(e(1)-en(1)^2*sin(theta)^2)/e(1); qn=sqrt(e(end)-en(1)^2*sin(theta)^2)/e(end); for j=2:(length(e)-1) beta=d(j)*2*pi/lambda*sqrt(e(j)-en(1)^2*sin(theta)^2); q=sqrt(e(j)-en(1)^2*sin(theta)^2)/e(j); em(j,1,1)=cos(beta); em(j,1,2)=-fukso*sin(beta)/q; em(j,2,1)=-fukso*sin(beta)*q; em(j,2,2)=cos(beta); end emtot=[1 0; 0 1]; for j=2:(length(e)-1) emtot1(:,:)=em(j,:,:); emtot=emtot*emtot1; end rp=((emtot(1,1)+emtot(1,2)*qn)*q1-(emtot(2,1)+emtot(2,2)*qn))/... ((emtot(1,1)+emtot(1,2)*qn)*q1+(emtot(2,1)+emtot(2,2)*qn)); tp=2*q1/((emtot(1,1)+emtot(1,2)*qn)*q1+(emtot(2,1)+emtot(2,2)*qn)); ref=rp*conj(rp); tra=tp*conj(tp)/cos(theta)*en(1)*qn; REF(jtheta)=ref; TRA(jtheta)=tra; end plot(THETA_ext_deg,REF,expdata(:,1),REF(1)/expdata(1,2)*expdata(:,2),'xk') ax=axis; text(ax(1)+(ax(2)-ax(1))*.1,ax(3)+(ax(4)-ax(3))*.2,['PRISM: n=',num2str(en(1)),... '\newline GOLD: n=',num2str(en(2)),'+i',num2str(ek(2)),' d=',num2str(d(2)*1e9),'nm',... '\newline LAYER: n=',num2str(en(3)),'+i',num2str(ek(3)),' d=',num2str(d(3)*1e9),'nm']); title(FNAME) pause(.05) end