Відкрити головне меню

Абендфарт

Приєднався 31 липня 2010
нема опису редагування
#include <stdio.h>
Абенфарт значить вечірня мандрівка
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_integration.h>
 
double sphbess(int n, double x)
{
double jb, jb1, t;
int i;
if (x>=2)
{
jb = sin(x)/x; jb1 = cos(x)/x;
for (i = 1; i <= n; i++)
{
t = jb;
jb = jb*(2*i-1)/x-jb1;
jb1 = t;
}
}
else
{
t = 1.0;
i = 0;
jb = 1;
jb1 = x*x/2;
while (jb != jb+t)
{
i++;
t *= -jb1/i/(2*n+2*i+1);
jb += t;
}
for (i=1; i <= n; i++) jb *= x/(2*i+1);
}
return jb;
}
 
double Lorential (double x, double Ga)
{
return Ga/M_PI/ (x*x + Ga*Ga);
}
 
typedef struct
{
int L;
double Q;
} istrc;
 
double f (double x, void * params)
{
double Q;
int n;
istrc *ad = (istrc *) params;
Q= ad->Q;
n = ad->L;
if (x <=1) return pow(x, n+2)*sphbess(n, Q*x);
if (x >1) return pow(x, -n+1)*sphbess(n, Q*x);
}
 
 
int main()
{
double w_p = 4.0; // plasma frequency
double eps_m0 = 1.0; //
double eps_ext = 1.0;
double wL; // plasmon for L
double hbar = 6.63e-34;
double m0 = 9.1e-31;
double eV = 1.6e-19;
double R = 20e-9;
double alpha, Q;
double Ei = 1.e4; // energy of incident particle
double k;
double beta=1.e-4;
double Pi=3.14;
double G;
double w=2.0;
int L;
double EL;
double Ga=1.e12;
double A = 5000;
int w_dependence = 1;
L=1;
gsl_integration_workspace * ws
= gsl_integration_workspace_alloc (1000);
double result, error;
gsl_function F;
F.function = &f;
if (w_dependence)
{
for (w = 1.e12; w < 1.e15; w+=1.e12)
{
double RP = 0;
for (L =1; L < 10; L++)
{
EL=sqrt(L*w_p/(L*eps_m0+(L+1)*eps_ext));
k = sqrt(2*m0* Ei*eV/(hbar*hbar));
wL=EL*eV/hbar;
istrc is;
Q = sqrt(k*k-m0*EL*eV/(hbar*hbar)-k*k*cos(beta)*sqrt(1-(2*m0*EL*eV)/(hbar*hbar*k*k)))*R;
is.L = L;
is.Q = Q;
F.params = &is;
gsl_integration_qags (&F, 0, A, 0, 1e-7, 1000,
ws, &result, &error);
G=-eV*(2L+1)*sqrt(2L+1)*2*Pi*sqrt(R*R*R*R*R*wL*wL/(L*eps_m0+(L+1)*eps_ext))*sqrt(hbar/(2*wL));
RP += 2*Pi*G*G*result*result/hbar*Lorential(w-wL,Ga)/hbar;
}
printf("%f %g \n", w, RP);
}
}
else
{
w = wL;
for (alpha = 0; alpha < 0.001; alpha +=0.000001)
{
double RR = 0;
for (L =1; L < 10; L++)
{
EL=sqrt(L*w_p/(L*eps_m0+(L+1)*eps_ext));
k = sqrt(2*m0* Ei*eV/(hbar*hbar));
wL=EL*eV/hbar;
istrc is;
Q = sqrt(k*k-m0*EL*eV/(hbar*hbar)-k*k*cos(alpha)*sqrt(1-(2*m0*EL*eV)/(hbar*hbar*k*k)))*R;
is.L = L;
is.Q = Q;
F.params = &is;
gsl_integration_qags (&F, 0, A, 0, 1e-7, 1000,
ws, &result, &error);
G=-eV*(2L+1)*sqrt(2L+1)*2*Pi*sqrt(R*R*R*R*R*wL*wL/(L*eps_m0+(L+1)*eps_ext))*sqrt(hbar/(2*wL));
RR += 2*Pi*G*G*result*result/hbar*Lorential(w-wL,Ga)/hbar;
}
printf("%f %g %g \n", alpha, Q, RR);
}
}
return 0;
}
437

редагувань