0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 function [Tous,Tousp,Tousn]=Calc_Tous(file,linlat,sigmas)
0029
0030 if nargin <1
0031 help Calc_Tous
0032 return;
0033 end
0034
0035
0036
0037
0038 try
0039 [dummy s ax bx mux etax etaxp az bz muz etaz etazp] = ...
0040 textread(linlat,'%s %f %f %f %f %f %f %f %f %f %f %f','headerlines',4);
0041 catch
0042 errmsg = lasterr;
0043 if strfind(errmsg, 'File not found.');
0044 error('Input file %s not found \n Abort \n',linlat);
0045 else
0046 error('Unknown error %s \n', errmsg)
0047 end
0048 end
0049
0050 ex = 3.7E-9;
0051 ez = ex*0.01;
0052 se = 1.016E-3;
0053 gx = (1+ax.^2)./bx;
0054
0055 sx = sqrt(bx*ex+(se*etax).^2);
0056 sxp= sqrt(gx*ex+(se*etaxp).^2);
0057 sz = sqrt(bz*ez+(se*etaz).^2);
0058
0059 [dummy sd acen lost nom] = textread(file,'%d %f %f %f %s','headerlines',3);
0060 Nb_pts=length(sd)/2;
0061 s=sd(1:Nb_pts);
0062 acen_p=acen(1:Nb_pts); acen_n=acen(Nb_pts+1:end);
0063
0064
0065 if nargin <3
0066 ss = 3.6E-3;
0067 fprintf('Par defaut la longueur du paquet vaut : sigma_s= %f m\n',ss)
0068 else
0069 ss =sigmas;
0070 end
0071
0072 V = 8*pi*sqrt(pi).*sx.*sz.*ss;
0073 E = 2.75;
0074 Gamma = 5382.57324;
0075 c = 2.99792458E8;
0076 re = 2.81794092E-15;
0077 I0 = 1.2E-3;
0078 Qe = 1.6022E-19;
0079 L = 354.097;
0080 N = I0/c*L/Qe;
0081
0082
0083
0084
0085
0086
0087 FF = sqrt(pi)*re*re*c*N*3600;
0088 Tinvp = zeros(1,Nb_pts);
0089 Tinvn = zeros(1,Nb_pts);
0090
0091 for i=1:Nb_pts,
0092 Tinvp(i) = (FF*Cxi((acen_p(i)/Gamma/sxp(i))^2)/sxp(i) ...
0093 /Gamma^3/acen_p(i)^2/V(i));
0094 Tinvn(i) = (FF*Cxi((acen_n(i)/Gamma/sxp(i))^2)/sxp(i) ...
0095 /Gamma^3/acen_n(i)^2/V(i));
0096 end
0097
0098
0099 Tinv = (Tinvp+Tinvn)/2;
0100
0101 Tp=1./Tinvp;
0102 Tn=1./Tinvn;
0103 T =1./Tinv;
0104
0105
0106 s1=[0 ; s]; s2=[s ; 0]; ds=s2-s1;
0107 Tousp=1/sum(1./Tp.*ds(1:end-1)')*s(end);
0108 Tousn=1/sum(1./Tn.*ds(1:end-1)')*s(end);
0109
0110 Tous=2/(1/Tousp+1/Tousn);
0111
0112
0113
0114
0115 figure(80)
0116 subplot(2,1,1)
0117 plot(s,Tp,'.-');
0118 grid on
0119 xlabel('s (m)')
0120 ylabel('Tp (h)')
0121
0122 subplot(2,1,2)
0123 plot(s,Tn,'.-');
0124 grid on
0125 xlabel('s (m)')
0126 ylabel('Tn (h)')
0127
0128
0129
0130
0131 files='/home/nadolski/matlab/soleil/structure';
0132 struc=dlmread(files);
0133
0134
0135 acen=acen*100;
0136
0137 figure(81)
0138 clf
0139 plot([sd(1:Nb_pts); 2*sd(Nb_pts)-sd(Nb_pts:-1:1)],[acen(1:Nb_pts) ; acen(Nb_pts:-1:1)],'b.-');
0140 hold on
0141 plot([sd(Nb_pts+1:end); 2*sd(end)-sd(end:-1:1+Nb_pts)],[acen(Nb_pts+1:end) ; acen(end:-1:Nb_pts+1)],'b.-');
0142 plot(struc(:,1),struc(:,2)/2,'k-')
0143 axis([0 88.5 -7 7])
0144 grid on
0145 title('solamor2')
0146 xlabel('s (m)')
0147 ylabel('Acceptance en energie \epsilon (%)')
0148