source: PSPA/parmelaPSPA/trunk/distrirhom.f @ 465

Last change on this file since 465 was 12, checked in by lemeur, 12 years ago

parmela pspa initial

File size: 1.8 KB
Line 
1        subroutine distrirhom(rmax,npoints,flag,nhom)
2c-----------------------------------------------------------------------
3c       Author Jerome Gonichon
4c       for x-y in input10 ran=3  (homogeneous distribution)
5c-----------------------------------------------------------------------
6c
7        include 'param_sz.h'
8        include 'constcom.h'
9c
10        common/hom/xhom(imaa),yhom(imaa)
11        dimension       xhom1(imaa),yhom1(imaa),rhom(imaa)
12c-----------------------------------------------------------------------
13        rmaxsq=rmax**2
14        delthom=2*rmax/(sqrt(float(npoints))-1)
15        nhom=0
16        x=-rmax
17        y=-rmax
18        do 10 i=1,npoints
19 20     continue
20        if(x.lt.rmax)then
21                x=x+delthom
22        else
23                if(y.lt.rmax)then
24                        y=y+delthom
25                        x=-rmax
26                else
27                goto 10
28                endif
29        endif
30        rsq=x**2+y**2
31        if(rsq.gt.rmaxsq)then
32        goto 20
33        endif
34        xhom(i)=x
35        yhom(i)=y
36        rhom(i)=sqrt(xhom(i)**2+yhom(i)**2)
37        write(50,*)rhom(i)
38        nhom=nhom+1
39 10     continue
40        if(flag.eq.0)return
41c---sorting rhom ascendant
42        do 12 j=2,nhom
43        a=rhom(j)
44        b=xhom(j)
45        c=yhom(j)
46                do 11 i=j-1,1,-1
47                if(rhom(i).le.a)goto 13
48                rhom(i+1)=rhom(i)
49                xhom(i+1)=xhom(i)
50                yhom(i+1)=yhom(i)
51 11             continue
52        i=0
53 13     rhom(i+1)=a
54        xhom(i+1)=b
55        yhom(i+1)=c
56 12     continue       
57        if(flag.eq.1)return
58        if(flag.eq.2)then
59        do 60 i=1,nhom
60        xhom1(i)=xhom(nhom+1-i)
61        yhom1(i)=yhom(nhom+1-i)
62 60     continue       
63        do 50 i=1,nhom
64        xhom(i)=xhom1(i)
65        yhom(i)=yhom1(i)
66        write(48,*)xhom(i),yhom(i),rhom(i)
67 50     continue
68        return
69        endif
70        inc=0
71        n=0
72        do 30 i=1,nhom/2+3
73        xhom1(i+inc)=xhom(i)
74        xhom1(i+inc+4)=xhom(nhom+1-i)
75        yhom1(i+inc)=yhom(i)
76        yhom1(i+inc+4)=yhom(nhom+1-i)
77        n=n+1
78        if(n.eq.4)then
79                inc=inc+4
80                n=0
81        endif
82 30     continue
83        do 40 i=1,nhom
84        xhom(i)=xhom1(i)
85        yhom(i)=yhom1(i)
86        write(49,*)xhom(i),yhom(i),sqrt(xhom(i)**2+yhom(i)**2)
87 40     continue
88        return
89        end
90c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
Note: See TracBrowser for help on using the repository browser.