00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 void electron_velocity(void)
00031 {
00032 int i, n, iv, el;
00033 double *vertex;
00034 double *c;
00035 double *aux;
00036 double *dum;
00037 double ksquared;
00038 double xvelocity;
00039 double yvelocity;
00040 double zvelocity;
00041 double thesquareroot;
00042 double superparticle_energy;
00043
00044 dum=malloc((Ng+1)*sizeof(double));
00045 vertex=malloc((Ne+1)*sizeof(double));
00046 aux=malloc((Ne+1)*sizeof(double));
00047 c=malloc((Ng+1)*sizeof(double));
00048 if(vertex==NULL || aux==NULL || c==NULL || dum==NULL){
00049 printf("electron_velocity : not enough memory in allocating arrays.\n");
00050 exit(0);
00051 }
00052
00053 for(i=0;i<Ne;i++){
00054 aux[i]=0.;
00055 vertex[i]=0.;
00056 }
00057
00058
00059
00060 for(n=1;n<=INUM;n++){
00061 iv=(int)(P[n][0]);
00062 el=(int)(P[n][9]);
00063 ksquared=P[n][1]*P[n][1]+P[n][2]*P[n][2]+P[n][3]*P[n][3];
00064 thesquareroot=sqrt(1.+4.*alphaK[i_dom[el]][iv]*HHM[i_dom[el]][iv]*ksquared);
00065 xvelocity=P[n][1]*HM[i_dom[el]][iv]/thesquareroot;
00066 vertex[el]+=1.;
00067 aux[el]+=xvelocity;
00068 }
00069
00070
00071
00072 for(i=0;i<Ne;i++) if(vertex[i]!=0.) aux[i]/=vertex[i];
00073 for(i=0;i<Ng;i++){ c[i]=0.; dum[i]=0.;}
00074 for(i=0;i<Ne;i++){
00075 int j;
00076 for(j=0;j<4;j++){
00077 dum[noeud_geo[j][i]-1]+=aux[i];
00078 c[noeud_geo[j][i]-1]+=1.;
00079 }
00080 }
00081 for(i=0;i<Ng;i++) dum[i]/=c[i];
00082 for(i=0;i<Ng;i++) XVEL[i]+=dum[i];
00083
00084
00085
00086 for(i=0;i<Ne;i++){
00087 aux[i]=0.;
00088 vertex[i]=0.;
00089 }
00090 for(n=1;n<=INUM;n++){
00091 iv=(int)(P[n][0]);
00092 el=(int)(P[n][9]);
00093 ksquared=P[n][1]*P[n][1]+P[n][2]*P[n][2]+P[n][3]*P[n][3];
00094 thesquareroot=sqrt(1.+4.*alphaK[i_dom[el]][iv]*HHM[i_dom[el]][iv]*ksquared);
00095 yvelocity=P[n][2]*HM[i_dom[el]][iv]/thesquareroot;
00096 vertex[el]+=1.;
00097 aux[el]+=yvelocity;
00098 }
00099
00100
00101
00102 for(i=0;i<Ne;i++) if(vertex[i]!=0.) aux[i]/=vertex[i];
00103 for(i=0;i<Ng;i++){ c[i]=0.; dum[i]=0.;}
00104 for(i=0;i<Ne;i++){
00105 int j;
00106 for(j=0;j<4;j++){
00107 dum[noeud_geo[j][i]-1]+=aux[i];
00108 c[noeud_geo[j][i]-1]+=1.;
00109 }
00110 }
00111 for(i=0;i<Ng;i++) dum[i]/=c[i];
00112 for(i=0;i<Ng;i++) YVEL[i]+=dum[i];
00113
00114
00115
00116
00117 for(i=0;i<Ne;i++){
00118 aux[i]=0.;
00119 vertex[i]=0.;
00120 }
00121 for(n=1;n<=INUM;n++){
00122 iv=(int)(P[n][0]);
00123 el=(int)(P[n][9]);
00124 ksquared=P[n][1]*P[n][1]+P[n][2]*P[n][2]+P[n][3]*P[n][3];
00125 thesquareroot=sqrt(1.+4.*alphaK[i_dom[el]][iv]*HHM[i_dom[el]][iv]*ksquared);
00126 zvelocity=P[n][3]*HM[i_dom[el]][iv]/thesquareroot;
00127 vertex[el]+=1.;
00128 aux[el]+=zvelocity;
00129 }
00130
00131
00132
00133 for(i=0;i<Ne;i++) if(vertex[i]!=0.) aux[i]/=vertex[i];
00134 for(i=0;i<Ng;i++){ c[i]=0.; dum[i]=0.;}
00135 for(i=0;i<Ne;i++){
00136 int j;
00137 for(j=0;j<4;j++){
00138 dum[noeud_geo[j][i]-1]+=aux[i];
00139 c[noeud_geo[j][i]-1]+=1.;
00140 }
00141 }
00142 for(i=0;i<Ng;i++) dum[i]/=c[i];
00143 for(i=0;i<Ng;i++) ZVEL[i]+=dum[i];
00144
00145
00146
00147 for(i=0;i<Ne;i++){
00148 aux[i]=0.;
00149 vertex[i]=0.;
00150 }
00151 for(n=1;n<=INUM;n++){
00152 iv=(int)(P[n][0]);
00153 el=(int)(P[n][9]);
00154 ksquared=P[n][1]*P[n][1]+P[n][2]*P[n][2]+P[n][3]*P[n][3];
00155 thesquareroot=sqrt(1.+4.*alphaK[i_dom[el]][iv]*HHM[i_dom[el]][iv]*ksquared);
00156 superparticle_energy=(thesquareroot-1.)/(2.*alphaK[i_dom[el]][iv]);
00157 if(iv==2) superparticle_energy+=EG[i_dom[el]];
00158 vertex[el]+=1.;
00159 aux[el]+=superparticle_energy;
00160 if(iv==2) aux[el] += EG[i_dom[el]];
00161 }
00162
00163
00164
00165 for(i=0;i<Ne;i++) if(vertex[i]!=0.) aux[i]/=vertex[i];
00166 for(i=0;i<Ng;i++){ c[i]=0.; dum[i]=0.;}
00167 for(i=0;i<Ne;i++){
00168 int j;
00169 for(j=0;j<4;j++){
00170 dum[noeud_geo[j][i]-1]+=aux[i];
00171 c[noeud_geo[j][i]-1]+=1.;
00172 }
00173 }
00174 for(i=0;i<Ng;i++) dum[i]/=c[i];
00175 for(i=0;i<Ng;i++) ENEL[i]+=dum[i];
00176
00177 free(c);
00178 free(aux);
00179 free(vertex);
00180 free(dum);
00181 }