1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
/*******************************************************************************
* *
* Viewmol *
* *
* I N S I N T . C *
* *
* Copyright (c) Joerg-R. Hill, December 2000 *
* *
********************************************************************************
*
* $Id: insint.c,v 1.5 2000/12/10 15:09:30 jrh Exp $
* $Log: insint.c,v $
* Revision 1.5 2000/12/10 15:09:30 jrh
* Release 2.3
*
* Revision 1.4 1999/05/24 01:26:03 jrh
* Release 2.2.1
*
* Revision 1.3 1999/02/07 21:50:39 jrh
* Release 2.2
*
* Revision 1.2 1998/01/26 00:48:09 jrh
* Release 2.1
*
* Revision 1.1 1996/12/10 18:41:30 jrh
* Initial revision
*
*/
#include<math.h>
#include "viewmol.h"
extern struct WINDOW windows[];
extern struct MOLECULE *molecules;
extern double temp;
void makeINSIntensity(void)
{
struct MOLECULE *mol;
double pi, h, c, bk, amu, cotfac, qfac;
register double f, rll, smax;
register int i, j, k;
/* This subroutine computes inelastic neutron scattering intensities */
mol=&molecules[windows[SPECTRUM].set];
pi=4.0*atan(1.0);
h=6.626176e-34;
c=2.99792458e10;
bk=1.380662e-23;
amu=1.6605655e-27;
if (temp > 0.0)
cotfac=h*c/(2.0*bk*temp);
else
cotfac=h*c/(2.0*bk*0.001);
qfac=h*1.e20/(4.0*pi*pi*c*amu);
smax=0.0;
for (i=0; i<3*mol->na; i++)
{
mol->normal_modes[i].ins_intensity=0.0;
if (mol->normal_modes[i].wavenumber > 0.0)
{
f=sqrt(qfac/(2.0*mol->normal_modes[i].wavenumber)
/tanh(cotfac*mol->normal_modes[i].wavenumber));
k=0;
for (j=0; j<mol->na; j++)
{
rll=sqrt(mol->cnm[i+mol->nmodes*k]*mol->cnm[i+mol->nmodes*k]
+mol->cnm[i+mol->nmodes*(k+1)]*mol->cnm[i+mol->nmodes*(k+1)]
+mol->cnm[i+mol->nmodes*(k+2)]*mol->cnm[i+mol->nmodes*(k+2)]);
k+=3;
mol->normal_modes[i].ins_intensity+=mol->atoms[j].neutronScatterfac*rll;
}
mol->normal_modes[i].ins_intensity*=f;
smax=smax > mol->normal_modes[i].ins_intensity
? smax : mol->normal_modes[i].ins_intensity;
}
}
for (i=0; i<3*mol->na; i++)
{
if (smax == 0.0)
mol->normal_modes[i].rel_ins_intensity=0.0;
else
mol->normal_modes[i].rel_ins_intensity=100.0*mol->normal_modes[i].ins_intensity/smax;
}
}
|
