Next: Appendice C
Up: Appendice B
Previous: Appendice B
Contents
Figure 7:
Diagramma di flusso del programma CURRECURRE
 |
c---------------------------------------------------------
c
c PROGRAM CURRECURRE
c
c---------------------------------------------------------
c INIZIALIZZAZIONE VARIABILI E PARAMETRI
parameter (nc=11)
c nc numero colonne del catalogo
parameter (na=236)
c na numero di oggetti del catalogo
parameter (q=1000000)
c q estrazioni casuali totali
parameter (p=726)
c p serve a cambiare il campione casuale
parameter (z=734231)
c z campione casuale al posto del catalogo
character*800 ivet
real a(na,nc)
real x(na,na),xx(na,na)
parameter (nr=100)
c PARAMETRI PER SELEZIONE OGGETTI CON
c PARTICOLARI PROPRIETA' FISICHE
c nr numero di intervalli
c parameter (jj=7)
c jj=4 magnitudini; jj=7 elongazioni
c parameter (kk=4)
c kk come jj
c parameter (xxx=100)
c xxx elongazione-magnitudine massima selezionata
c parameter (yyy=24.299)
c yyy elongazione-magnitudine minima selezionata
real xo
parameter (qw=15000)
c qw numero di oggetti del campione casuale
real zzz(qw,3) ,zz(q,3) ,bv(nr,7),ls(nr,2)
real xw(qw,qw)
real aa(nr) ,suma
real ab(nr) ,sumb
real ac(nr) ,sumc
real xk(na,qw)
c AZZERAMENTO VARIABILI E MATRICI
do i=1,qw
do j=1,qw
xw(i,j)=0
enddo
do j=1,3
zzz(i,j)=0
zz(i,j)=0
enddo
do j=1,na
xk(j,i)=0
enddo
enddo
do i=1,q
zz(i,1)=0
zz(i,2)=0
zz(i,3)=0
enddo
do i=1,nr
aa(i)=0
ab(i)=0
ac(i)=0
ls(i,1)=0
ls(i,2)=0
enddo
do i=1,na
do j =1,nc
a(i,j) = 0
enddo
do k =1,na
x(i,k) = 0
xx(i,k) = 0
enddo
enddo
c ACQUISIZIONE CATALOGO
open(16,file='Cat400I.dat',status='old')
do 10 k= 1,16
read(16,'(a800)',end=2) ivet
10 continue
read(16,*,end=2) ((a(i,j),j=1,nc),i=1,na)
2 close(16)
c SELEZIONE OGGETTI PER PROPRIETA' FISICHE PARTICOLARI
c
c write (*,*) a(1,1) ,a(1,2) ,a(2,1)
c do 999 i=1,na
c if (a(i,kk).lt.yyy.and.a(i,jj).lt.xxx) then
c a(i,2)=a(i,2)
c a(i,3)=a(i,3)
c else
c a(i,2)=0.0
c a(i,3)=0.0
c endif
c999 continue
c CALCOLO ESTREMI CAMPO DI SIMULAZIONE
xmin=10000000
xmax=-100000000
ymin=10000000
ymax=-100000000
zmin=10000000
zmax=-100000000
do 234 i=1,na
if (a(i,2).le.xmin) then
xmin = a(i,2)
elseif (a(i,2).ge.xmax) then
xmax = a(i,2)
endif
234 continue
do 235 i=1,na
if (a(i,3).le.ymin) then
ymin = a(i,3)
elseif (a(i,3).ge.ymax) then
ymax = a(i,3)
endif
235 continue
do 238 i=1,na
if (a(i,4).le.zmin) then
zmin = a(i,4)
elseif (a(i,4).ge.zmax) then
zmax = a(i,4)
endif
238 continue
write(*,300) xmin ,xmax ,ymin ,ymax ,zmin ,zmax
c Estrazione campione casuale di
c oggetti uniformemente distribuiti
do 20001 i = 1,q
do 22001 k =1,3
zz(i,k)=ran3(idum)
22001 continue
20001 continue
c CICLO PER SOSTITUZIONE CATALOGO CON CAMPIONE CASUALE
c PER MISURA ERRORI SISTEMATICI DELLO STIMATORE
c
c do 21001 i=1,na
c
c
c if (a(i,2).ne.0.0) then
c
c a(i,2)=(xmax - xmin)*zz(i+z,1)+ xmin
c a(i,3)=(ymax - ymin)*zz(i+z,2)+ ymin
c
c endif
c
c
c21001 continue
c open(80,file='celongcasu.dat',status='unknown')
c write(80,127) ((a(i,j),j=1,nc),i=1,na)
c127 format(15g20.10)
c close(80)
c ------------- CICLO DI CALCOLO ----------
do 20 i=1,na
do 30 j =1,na
c if (a(i,2).ne.0.and.a(j,2).ne.0) then
ax=((a(i,2)-a(j,2))**2)
ay=((a(i,3)-a(j,3))**2)
az=((a(i,4)-a(j,4))**2)
x(i,j)=sqrt(ax+ay+az)
c else
c x(i,j)=0
c endif
30 continue
20 continue
c do 21 i=1,na
c do 31 j =1,na
c
c if (xx(i,j).ne.0) then
c x(i,j)=log(xx(i,j))/2.3022585093
c else
c x(i,j)=0
c endif
c
c31 continue
c21 continue
rmin=10000000
rmax=-100000000
do 236 i=1,na
do 237 j=1,na
if (x(i,j).le.rmin.and.x(i,j).ne.0) then
rmin = x(i,j)
elseif (x(i,j).ge.rmax) then
rmax = x(i,j)
endif
237 continue
236 continue
c rmax= 8000.0
c rmin= 0.0
write(*,300) rmin ,rmax
xo = (rmax -rmin)/nr
write(*,300) xo
do 23 i=1,na
do 34 j=1,na
do 37 n=1,nr
if (x(i,j).le.rmin+n*xo.and.x(i,j).gt.rmin+(n-1)*xo)then
aa(n) = aa(n) + 1
endif
37 continue
34 continue
23 continue
aa(1) = aa(1) + 2
suma = 0
do 47 i = 1,nr
write(*,300) aa(i)
suma = suma + aa(i)
47 continue
write(*,300) suma
c do 2000 i = 1,q
c do 2200 k =1,2
c
c zz(i,k)=ran3(idum)
c
c2200 continue
c2000 continue
do 2100 i=1,qw
zzz(i,1)=(xmax - xmin)*zz(i+p,1)+ xmin
zzz(i,2)=(ymax - ymin)*zz(i+p,2)+ ymin
zzz(i,3)=(zmax - zmin)*zz(i+p,3)+ zmin
2100 continue
do 2002 i=1,qw
do 3002 j =1,qw
a1=((zzz(i,1)-zzz(j,1))**2)
a2=((zzz(i,2)-zzz(j,2))**2)
a3=((zzz(i,3)-zzz(j,3))**2)
c if(sqrt(a1+a2+a3).ne.0)then
xw(i,j)=sqrt(a1+a2+a3)
c xw(i,j)=log(sqrt(a1+a2))/2.3022585093
c else
c xw(i,j)=0
c endif
3002 continue
2002 continue
do 231 i=1,qw
do 341 j=1,qw
do 371 n=1,nr
if (xw(i,j).le.rmin+n*xo.and.xw(i,j).gt.rmin+(n-1)*xo) then
ab(n) = ab(n) + 1
endif
371 continue
341 continue
231 continue
sumb = 0
do 471 i = 1,nr
write(*,300) ab(i)
sumb = sumb + ab(i)
471 continue
write(*,300) sumb
do 203 i=1,na
do 303 j =1,qw
a4=((a(i,2)-zzz(j,1))**2)
a5=((a(i,3)-zzz(j,2))**2)
a6=((a(i,4)-zzz(j,3))**2)
c if(a(i,2).ne.0.and.sqrt(a3+a4).ne.0) then
xk(i,j)=sqrt(a4+a5+a6)
c xk(i,j)=log(sqrt(a3+a4))/2.3022585093
c else
c xk(i,j)=0
c endif
303 continue
203 continue
do 2361 i=1,na
do 346 j=1,qw
do 376 n=1,nr
if (xk(i,j).le.rmin+n*xo.and.xk(i,j).gt.rmin+(n-1)*xo) then
ac(n) = ac(n) + 1
endif
376 continue
346 continue
2361 continue
sum = 0
do 478 i = 1,nr
write(*,300) ac(i)
sumc = sumc + ac(i)
478 continue
write(*,300) sumc
c --------- CALCOLO STIMATORI -----------
do 13i=1,nr
bv(i,1)=aa(i)/(suma)
bv(i,2)=ab(i)/(sumb)
bv(i,3)=rmin+i*xo
bv(i,4)=bv(i,1)/bv(i,2) -1
bv(i,5)=ac(i)/(sumc)
bv(i,6)=(bv(i,1)+bv(i,2)-2*bv(i,5))/bv(i,2)
bv(i,7)=bv(i,1)/bv(i,5) -1
ls(i,1)=bv(i,3)
ls(i,2)=bv(i,6)
13 continue
c ---------- FILE DI OUTPUT -----------
c FILE CON DIVERSI STIMATORI PER VISUALIZZAZIONE SM
open(404,file='Cat400I.cor',status='unknown')
write(404,1253) ((bv(i,j),j=1,7),i=1,nr)
1253 format(7g20.10)
close(404)
c FILE CON STIMATORE LS PER STIMA PARAMETRI DA FIT
open(4043,file='Cat400Ils',status='unknown')
write(4043,1223) ((ls(i,j),j=1,2),i=1,nr)
1223 format(2g20.10)
close(4043)
300 format(6g16.6)
301 format(10i10)
stop
end
c ROUTINE PER ESTRAZIONE CAMPIONI CASUALI
c DA NUMERICAL RECIPES
function ran2(idum)
integer idum,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
real ran2,am,eps,rnmx
parameter (im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1)
parameter (ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211)
parameter (ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7)
parameter (rnmx=1.-eps)
integer idum2,j,k,iv(ntab),iy
save iv,iy,idum2
data idum2/123456789/, iv/ntab*0/, iy/0/
if (idum.le.0) then
idum=max(-idum,1)
idum2=idum
do j=ntab+8,1,-1
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
if (idum.lt.0) idum=idum+im1
if (j.le.ntab) iv(j)=idum
enddo
iy=iv(1)
endif
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
if(idum.lt.0) idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
if(idum2.lt.0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1) iy=iy+imm1
ran2=min(am*iy,rnmx)
return
end
function ran3(idum)
integer idum
real mbig,mseed,mz
real ran3,fac
parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=1./mbig)
integer i,iff,ii,inext,inextp,k
real mj,mk,ma(55)
save iff,inext,inextp,ma
data iff /0/
if (idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz)mk=mk+mbig
mj=ma(ii)
enddo
do k=1,4
do i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
enddo
enddo
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran3=mj*fac
return
end
Maurilio Pannella
2001-07-30