Line | |
---|
1 | subroutine distriphase(xmax,sigma,npoints) |
---|
2 | c for phase in input10 ran=2 and 3 (LAL) |
---|
3 | c----------------------------------------------------------------------- |
---|
4 | c |
---|
5 | include 'param_sz.h' |
---|
6 | include 'constcom.h' |
---|
7 | c |
---|
8 | common/dist/distp(imaa),xr(imaa),yr(imaa) |
---|
9 | dimension n(imaa/5) |
---|
10 | c----------------------------------------------------------------------- |
---|
11 | c* |
---|
12 | sqtwopi=sqrt(twopi) |
---|
13 | fac=xmax/sigma |
---|
14 | if(fac.lt.0.5)then |
---|
15 | do 10, i=1,npoints |
---|
16 | distp(i)=-xmax+(i-1)*2*xmax/npoints |
---|
17 | distp(i)=-distp(i) |
---|
18 | 10 continue |
---|
19 | else |
---|
20 | ech=xmax/1.4142/sigma |
---|
21 | const=npoints/sqtwopi/erf(ech)/sigma |
---|
22 | sigsq=2*sigma**2 |
---|
23 | nbin=int(min(npoints/5,100)) |
---|
24 | if (mod(nbin,2).ne.0.)then |
---|
25 | nbin=nbin+1 |
---|
26 | endif |
---|
27 | nhbin=nbin/2 |
---|
28 | pas=xmax/nhbin |
---|
29 | k=0 |
---|
30 | ntot=0 |
---|
31 | do 20,i=1,nhbin |
---|
32 | n(i)=int(const*exp(-((nhbin-i+0.5)*pas)**2/sigsq)*pas+0.5) |
---|
33 | ntot=ntot+n(i) |
---|
34 | 20 continue |
---|
35 | if(ntot.lt.int(npoints/2)+1)then |
---|
36 | nextra=int(npoints/2)+1-ntot |
---|
37 | do 25,l=nhbin,nhbin-nextra+1,-1 |
---|
38 | n(l)=n(l)+1 |
---|
39 | 25 continue |
---|
40 | endif |
---|
41 | do 30,i=1,nhbin |
---|
42 | do 35,j=k+1,k+n(i) |
---|
43 | if(j.gt.npoints/2+1)then |
---|
44 | nj=j |
---|
45 | go to 40 |
---|
46 | endif |
---|
47 | distp(j)=(n(i)-j+k)*pas/n(i)+(nhbin-i)*pas |
---|
48 | 35 continue |
---|
49 | k=k+n(i) |
---|
50 | 30 continue |
---|
51 | 40 continue |
---|
52 | if(nj.eq.0) nj=j |
---|
53 | do 50, i=nj,npoints |
---|
54 | distp(i)=-distp(2*nj-i-2) |
---|
55 | 50 continue |
---|
56 | endif |
---|
57 | return |
---|
58 | end |
---|
59 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|
Note: See
TracBrowser
for help on using the repository browser.