test_plasmadisp_jd

 test_plasmadisp_jd

0001 % test_plasmadisp_jd
0002 %
0003 clear all
0004 close all
0005 %
0006 p_opt = -1;
0007 %
0008 nx = 101;
0009 nx_as = 21;
0010 x = linspace(-10,10,nx);
0011 x_as = logspace(-10,10,nx_as);
0012 NZ_list = [10,10,100,1000,10000];
0013 NZn = length(NZ_list);
0014 %
0015 % Matlab function
0016 %
0017 Z_mat = NaN*ones(NZn,nx);
0018 Zp_mat = NaN*ones(NZn,nx);
0019 Zpp_mat = NaN*ones(NZn,nx);
0020 time_mat = NaN*ones(1,NZn);
0021 %
0022 Z_as_mat = NaN*ones(NZn,nx_as);
0023 Zp_as_mat = NaN*ones(NZn,nx_as);
0024 Zpp_as_mat = NaN*ones(NZn,nx_as);
0025 %
0026 leg0 = {};
0027 %
0028 for iZ = 1:NZn,
0029     tic
0030     [Z_mat(iZ,:),Zp_mat(iZ,:),Zpp_mat(iZ,:)] = Zplasma(x,NZ_list(iZ));
0031     time_mat(iZ) = toc;
0032     [Z_as_mat(iZ,:),Zp_as_mat(iZ,:),Zpp_as_mat(iZ,:)] = Zplasma(x_as,NZ_list(iZ));
0033     leg0{iZ} = ['mat, N = ',num2str(NZ_list(iZ))];
0034 end
0035 %
0036 leg0 = leg0(2:end);
0037 %
0038 % erfcc
0039 %
0040 tic
0041 xx = [x,zeros(1,1024 - nx)];
0042 %
0043 [yr,yi] = erfcc(nx,real(xx),imag(xx));
0044 %
0045 Z_erf = sqrt(pi)*(i*yr(1:nx) - yi(1:nx));
0046 Zp_erf = -2*(1 + x.*Z_erf);
0047 Zpp_erf = -2*(Z_erf - 2*x.*(1 + x.*Z_erf));
0048 %
0049 time_erf = toc;
0050 %
0051 xx_as = [x_as,zeros(1,1024 - nx_as)];
0052 %
0053 [yr_as,yi_as] = erfcc(nx_as,real(xx_as),imag(xx_as));
0054 %
0055 Z_as_erf = sqrt(pi)*(i*yr_as(1:nx_as) - yi_as(1:nx_as));
0056 Zp_as_erf = -2*(1 + x_as.*Z_as_erf);
0057 Zpp_as_erf = -2*(Z_as_erf - 2*x_as.*(1 + x_as.*Z_as_erf));
0058 %
0059 % series
0060 %
0061 Z_as_ser = NaN*ones(1,nx_as);
0062 Zp_as_ser = NaN*ones(1,nx_as);
0063 Zpp_as_ser = NaN*ones(1,nx_as);
0064 theory = NaN*ones(1,nx_as);
0065 theoryp = NaN*ones(1,nx_as);
0066 theorypp = NaN*ones(1,nx_as);
0067 %
0068 nx1 = (nx_as + 1)/2;
0069 %
0070 xpow_as = x_as(1:nx1-1);
0071 %
0072 Z_as_ser(1:nx1-1) = i*sqrt(pi)*exp(-xpow_as.^2) - 2*xpow_as.*...
0073     (1 - 2*xpow_as.^2/3 + 4*xpow_as.^4/15);% - 8*xpow_as.^6/105
0074 Zp_as_ser(1:nx1-1) = -2*i*xpow_as*sqrt(pi).*exp(-xpow_as.^2) - 2.*...
0075     (1 - 2*xpow_as.^2 + 4*xpow_as.^4/3);% - 8*xpow_as.^6/15
0076 Zpp_as_ser(1:nx1-1) = 2*i*sqrt(pi)*(2*xpow_as.^2 - 1).*exp(-xpow_as.^2) + 8*xpow_as.*...
0077     (1 - 4*xpow_as.^2/3);% + 4*xpow_as.^4/5
0078 %
0079 theory(1:nx1-1) = (16*xpow_as.^7/105)./(2*xpow_as);
0080 theoryp(1:nx1-1) = (16*xpow_as.^6/15)./(2);
0081 theorypp(1:nx1-1) = (32*xpow_as.^5/5)./(8*xpow_as);
0082 %
0083 theory(1:nx1-1) = max([theory(1:nx1-1);eps*ones(1,nx1-1)]);
0084 theoryp(1:nx1-1) = max([theoryp(1:nx1-1);eps*ones(1,nx1-1)]);
0085 theorypp(1:nx1-1) = max([theorypp(1:nx1-1);eps*ones(1,nx1-1)]);
0086 %
0087 xinv_as = 1./x_as(nx1+1:end);
0088 %
0089 Z_as_ser(nx1+1:end) = i*sqrt(pi)*exp(-1./xinv_as.^2) - xinv_as.*...
0090     (1 + xinv_as.^2/2 + 3*xinv_as.^4/4);% + 15*xinv_as.^6/8
0091 Zp_as_ser(nx1+1:end) = -2*i./xinv_as*sqrt(pi).*exp(-1./xinv_as.^2) + xinv_as.^2.*...
0092     (1 + 3*xinv_as.^2/2 + 15*xinv_as.^4/4);% + 105*xinv_as.^6/8
0093 Zpp_as_ser(nx1+1:end) = 2*i*sqrt(pi)*(2./xinv_as.^2 - 1).*exp(-1./xinv_as.^2) - 2*xinv_as.^3.*...
0094     (1 + 3*xinv_as.^2 + 45*xinv_as.^4/4);% + 105*xinv_as.^6/2
0095 %
0096 theory(nx1+1:end) = (15*xinv_as.^7/8)./(xinv_as);
0097 theoryp(nx1+1:end) = (105*xinv_as.^8/8)./(xinv_as.^2);
0098 theorypp(nx1+1:end) = (105*xinv_as.^9)./(2*xinv_as.^3);
0099 %
0100 theory(nx1+1:end) = max([theory(nx1+1:end);eps*ones(1,nx1-1)]);
0101 theoryp(nx1+1:end) = max([theoryp(nx1+1:end);2*eps./xinv_as.^2]);
0102 theorypp(nx1+1:end) = max([theorypp(nx1+1:end);4*eps./xinv_as.^4]);
0103 %
0104 Z_thei = sqrt(pi)*exp(-x.^2);
0105 Zp_thei = -2*x*sqrt(pi).*exp(-x.^2);
0106 Zpp_thei = 2*sqrt(pi)*(2*x.^2 - 1).*exp(-x.^2);
0107 %
0108 % figures
0109 %
0110 figure(1),clf
0111 %
0112 xlim = [-10,10];
0113 ylim = NaN;
0114 %
0115 leg = ['erfcc',leg0,'theory'];
0116 %
0117 graph1D_jd(x,real(Z_erf),0,0,'x','Z(x)','',NaN,xlim,ylim,'-','none','k',2,20,gca,0.9,0.7,0.7);
0118 graph1D_jd(x,real(Z_mat(2:end,:)).',0,0,'','','',NaN,xlim,ylim,'-','none',NaN,2);
0119 graph1D_jd(x,Z_thei,0,0,'','','',leg,xlim,ylim,'-','none','m',2);
0120 graph1D_jd(x,imag(Z_erf),0,0,'','','',NaN,xlim,ylim,'--','none','k',2);
0121 graph1D_jd(x,imag(Z_mat(2:end,:)).',0,0,'','','',NaN,xlim,ylim,'-','none',NaN,0.5);
0122 %
0123 figure(2),clf
0124 %
0125 xlim = [1,100000];
0126 ylim = NaN;
0127 %
0128 leg = {'erfcc','mat'};
0129 %
0130 graph1D_jd(xlim,time_erf*[1,1],1,1,'N','time (s)','',NaN,xlim,ylim,'--','none','k',2,20,gca,0.9,0.7,0.7);
0131 graph1D_jd(NZ_list(2:end),time_mat(2:end),1,1,'','','',leg,xlim,ylim,'none','+','r',2);
0132 %
0133 figure(3),clf,
0134 %
0135 xlim = [1e-10,1e10];
0136 ylim = NaN;%[1e-30,1];
0137 %
0138 leg = ['erfcc',leg0,'theory'];
0139 %
0140 graph1D_jd(x_as,abs(real(Z_as_erf) - real(Z_as_ser))./abs(real(Z_as_ser)),1,1,'x','Z(x)','',NaN,xlim,ylim,'none','o','k',2,20,gca,0.9,0.7,0.7);
0141 graph1D_jd(x_as,(abs(real(Z_as_mat(2:end,:)) - repmat(real(Z_as_ser),[NZn-1,1]))./abs(repmat(real(Z_as_ser),[NZn-1,1]))).',1,1,'','','',NaN,xlim,ylim,'none','+',NaN,2);
0142 graph1D_jd(x_as,theory,1,1,'','','',leg,xlim,ylim,'--','none',NaN,2);
0143 %
0144 figure(4),clf
0145 %
0146 xlim = [-10,10];
0147 ylim = NaN;
0148 %
0149 leg = ['erfcc',leg0,'theory'];
0150 %
0151 graph1D_jd(x,real(Zp_erf),0,0,'x','Z''(x)','',NaN,xlim,ylim,'-','none','k',2,20,gca,0.9,0.7,0.7);
0152 graph1D_jd(x,real(Zp_mat(2:end,:)).',0,0,'','','',NaN,xlim,ylim,'-','none',NaN,2);
0153 graph1D_jd(x,Zp_thei,0,0,'','','',leg,xlim,ylim,'-','none','m',2);
0154 graph1D_jd(x,imag(Zp_erf),0,0,'','','',NaN,xlim,ylim,'--','none','k',2);
0155 graph1D_jd(x,imag(Zp_mat(2:end,:)).',0,0,'','','',NaN,xlim,ylim,'-','none',NaN,0.5);
0156 %
0157 figure(5),clf,
0158 %
0159 xlim = [1e-10,1e10];
0160 ylim = NaN;%[1e-30,1];
0161 %
0162 leg = ['erfcc',leg0,'theory'];
0163 %
0164 graph1D_jd(x_as,abs(real(Zp_as_erf) - real(Zp_as_ser))./abs(real(Zp_as_ser)),1,1,'x','Z''(x)','',NaN,xlim,ylim,'none','o','k',2,20,gca,0.9,0.7,0.7);
0165 graph1D_jd(x_as,(abs(real(Zp_as_mat(2:end,:)) - repmat(real(Zp_as_ser),[NZn-1,1]))./abs(repmat(real(Zp_as_ser),[NZn-1,1]))).',1,1,'','','',NaN,xlim,ylim,'none','+',NaN,2);
0166 graph1D_jd(x_as,theoryp,1,1,'','','',leg,xlim,ylim,'--','none',NaN,2);
0167 %
0168 figure(6),clf
0169 %
0170 xlim = [-10,10];
0171 ylim = NaN;
0172 %
0173 leg = ['erfcc',leg0,'theory'];
0174 %
0175 graph1D_jd(x,real(Zpp_erf),0,0,'x','Z''''(x)','',NaN,xlim,ylim,'-','none','k',2,20,gca,0.9,0.7,0.7);
0176 graph1D_jd(x,real(Zpp_mat(2:end,:)).',0,0,'','','',NaN,xlim,ylim,'-','none',NaN,2);
0177 graph1D_jd(x,Zpp_thei,0,0,'','','',leg,xlim,ylim,'-','none','m',2);
0178 graph1D_jd(x,imag(Zpp_erf),0,0,'','','',NaN,xlim,ylim,'--','none','k',2);
0179 graph1D_jd(x,imag(Zpp_mat(2:end,:)).',0,0,'','','',NaN,xlim,ylim,'-','none',NaN,0.5);
0180 %
0181 figure(7),clf,
0182 %
0183 xlim = [1e-10,1e10];
0184 ylim = NaN;%[1e-30,1];
0185 %
0186 leg = ['erfcc',leg0,'theory'];
0187 %
0188 graph1D_jd(x_as,abs(real(Zpp_as_erf) - real(Zpp_as_ser))./abs(real(Zpp_as_ser)),1,1,'x','Z''''(x)','',NaN,xlim,ylim,'none','o','k',2,20,gca,0.9,0.7,0.7);
0189 graph1D_jd(x_as,(abs(real(Zpp_as_mat(2:end,:)) - repmat(real(Zpp_as_ser),[NZn-1,1]))./abs(repmat(real(Zpp_as_ser),[NZn-1,1]))).',1,1,'','','',NaN,xlim,ylim,'none','+',NaN,2);
0190 graph1D_jd(x_as,theorypp,1,1,'','','',leg,xlim,ylim,'--','none',NaN,2);
0191 %
0192 print_jd(p_opt,'Fig_Z',NaN,1);
0193 print_jd(p_opt,'Fig_t',NaN,2);
0194 print_jd(p_opt,'Fig_Zrel',NaN,3);
0195 print_jd(p_opt,'Fig_Zp',NaN,4);
0196 print_jd(p_opt,'Fig_Zprel',NaN,5);
0197 print_jd(p_opt,'Fig_Zpp',NaN,6);
0198 print_jd(p_opt,'Fig_Zpprel',NaN,7);
0199 %
0200 
0201 
0202