program ErrorCalculus
use moduleImageCreator
use MODEXCHG
use modnetcdf
implicit none
character(256)::parameterFileName,netcdfFileName,buff,imgN1,algo
integer::np,i,j,k,index,ib,intensity
real::x,y,dx,dy
real::h,d,r,lim,g,omega
integer::i0,i1,j0,j1
logical::civ1_,civ2_,patch1_,patch2_
real,dimension(:,:),allocatable::vect
integer::vec_size
real::meanX,meanY,rmsX,rmsY
integer::nbins
real,dimension(:,:),allocatable::errorX,errorY
real::decPart
integer::z,cptFalse
logical::convectFlow
integer::ik,jk
real::dx_tmp,dy_tmp
logical::display
real::minX,minY
minx=1.
miny=1.
civ1_=.false.
civ2_=.false.
patch1_=.false.
patch2_=.false.
convectFlow=.false.
nbins=10
np=iargc()
if(np==0) then
print*,"Usage: errorCalculus -f paramFile -p netcdfFile civ1-2 cv(if convectFlow)"
stop
endif
do i=1,np
CALL GETARG(i,buff)
if(buff=='-f') call GETARG(i+1,parameterFileName)
if(buff=='-p') call GETARG(i+1,netcdfFileName)
if(buff=='civ1') then
civ1_=.true.
algo=buff
endif
if(buff=='civ2') then
civ2_=.true.
algo=buff
endif
if(buff=='patch1') then
patch1_=.true.
algo=buff
endif
if(buff=='patch2') then
patch2_=.true.
algo=buff
endif
if(buff=='cv') then
convectFlow=.true.
endif
enddo
call readParametre_xml(parameterFileName)
BurgerCentreX=BurgerCentreX-1
print*,"start reading netcdf file",trim(adjustl(netcdfFileName))
call rd_netcdf(netcdfFileName)
if(civ1_)then
allocate(vect(nb_vectors,12))
vec_size=nb_vectors
vect(:,1:4)=V(:,1:4)
vect(:,9)=V(:,5)
elseif(patch1_) then
allocate(vect(nb_vec_patch,12))
vec_size=nb_vec_patch
vect(:,1:4)=VD(:,1:4)
0
elseif(civ2_) then
allocate(vect(nb_vectors2,12))
vec_size=nb_vectors2
vect(:,1:4)=V2(:,1:4)
vect(:,9)=V2(:,5)
elseif(patch2_) then
allocate(vect(nb_vec2_patch,12))
vec_size=nb_vec2_patch
vect(:,1:4)=VD2(:,1:4)
vect(:,9)=0
endif
vect(:,11)=vect(:,1)
vect(:,12)=vect(:,2)
!do ik=-100,100
!do jk=-100,100
!vect(:,1)=vect(:,11)+real(ik/10.)
!vect(:,2)=vect(:,12)+real(jk/10.)
if(convectFlow) then
print*,"Correcting Convect Flow"
vect(:,1)=vect(:,1)-vect(:,3)/2.
vect(:,2)=vect(:,2)-vect(:,4)/2.
endif
cptFalse=0
vect(:,10)=1
do i=1,vec_size
call burger(BurgerRadius,BurgerCentreX,BurgerCentreY,BurgerCirculation,dx,dy,vect(i,1),vect(i,2),omega)
vect(i,5)=dx+GlobalShiftX
vect(i,6)=dy+GlobalShiftY
vect(i,7)=vect(i,5)-vect(i,3)
vect(i,8)=vect(i,6)-vect(i,4)
! if((abs(vect(i,7))>(abs(vect(i,5)/10.))).or.(abs(vect(i,8))>(abs(vect(i,6)/10.)))) then
! print*,vect(i,5),vect(i,3),vect(i,6),vect(i,4)
! vect(i,7:10)=0.
! cptFalse=cptFalse+1
! endif
enddo
!if(cptFalse>150) goto 111
rmsX=0.
rmsY=0.
meanX=sum((vect(:,7)))/float(vec_size-cptFalse)
meanY=sum((vect(:,8)))/float(vec_size-cptFalse)
display=.false.
if(abs(meanX)<minx) then
minx=abs(meanX)
display=.true.
endif
if(abs(meanY)<miny) then
miny=abs(meanY)
display=.true.
endif
!if(display) then
do i=1,vec_size
if(vect(i,10)==1) then
rmsX=rmsX+sqrt( (vect(i,3)-vect(i,5))**2)
rmsY=rmsY+sqrt( (vect(i,4)-vect(i,6))**2)
endif
enddo
rmsX=rmsX/float(vec_size-cptFalse)
rmsY=rmsY/float(vec_size-cptFalse)
print*,"meanX",meanX
print*,"meanY",meanY
print*,"rmsX",rmsX
print*,"rmsY",rmsY
print*,"Mean Correlation",sum(vect(:,9))/float(vec_size-cptFalse)
print*,"Nb of False vectors: ",cptFalse," over ",vec_size,"..",real(ik/4.)," ",real(jk/4.)
!endif
111 continue
!enddo
!enddo
allocate(errorX(nbins,0:10),errorY(nbins,0:10))!,rmsX(nbins,0:10),rmsY(nbins,0:10),trueX(nbins,0:10),trueY(nbins,0:10))
errorX=0.
errorY=0.
!rmsX=0.
!rmsY=0.
!trueX=0.
!trueY=0.
! We put each vectors in a bin
do i=1,vec_size
if(vect(i,10)==1) then
decPart=abs(vect(i,5)-int(vect(i,5))) ! 0->1
z=1+int(decPart*(nbins)) !1->nbins
errorX(z,1)=errorX(z,1)+vect(i,7)
errorX(z,0)=errorX(z,0)+1
! print*,"vect(i,5) ",vect(i,5),"decPart ",decPart," bins: ",z
decPart=abs(vect(i,6)-int(vect(i,6))) ! 0->1
z=1+int(decPart*(nbins)) !1->nbins
errorY(z,1)=errorY(z,1)+vect(i,8)
errorY(z,0)=errorY(z,0)+1
endif
enddo
errorX(:,1)=errorX(:,1)/errorX(z,0)
errorY(:,1)=errorY(:,1)/errorY(z,0)
do i=1,vec_size
if(vect(i,10)==1) then
decPart=abs(vect(i,5)-int(vect(i,5))) ! 0->1
z=1+int(decPart*(nbins)) !1->nbins
errorX(z,2)=errorX(z,2)+(vect(i,7)-errorX(z,1))**2
errorX(z,3)=errorX(z,3)+(vect(i,7))**2
decPart=abs(vect(i,6)-int(vect(i,6))) ! 0->1
z=1+int(decPart*(nbins)) !1->nbins
errorY(z,2)=errorY(z,2)+(vect(i,8)-errorY(z,1))**2
errorY(z,3)=errorY(z,3)+(vect(i,8))**2
endif
enddo
open(unit=9,file=trim(adjustl(netcdfFileName))//'_Errors_'//trim(adjustl(algo))//'.dat',form='formatted',status='unknown')
do i=1,nbins
write(9,*) i/float(nbins)-1/(2.*nbins),errorX(i,1)/errorX(i,0),errorY(i,1)/errorY(i,0),sqrt(errorX(i,2)/errorX(i,0)),&
sqrt(errorY(i,2)/errorY(i,0)),sqrt(errorX(i,3)/errorX(i,0)),sqrt(errorY(i,3)/errorY(i,0))
write(*,*) i/float(nbins)-1/(2.*nbins),errorX(i,1)/errorX(i,0),errorY(i,1)/errorY(i,0),sqrt(errorX(i,2)/errorX(i,0)),&
sqrt(errorY(i,2)/errorX(i,0)),sqrt(errorX(i,3)/errorX(i,0)),sqrt(errorY(i,3)/errorX(i,0))
enddo
close(9)
end
!*--------------calculates the final position of a point subject
!*--------------to burgers vortex center at ccol,crow size r0 scale circ
!* Burgers' vortex on a regular grid
subroutine burger(r0,ccol,crow,circ,ddx,ddy,sx,sy,omeg)
real:: r0,xj,yj,rad,rad2,r02,a1,a2,a12,omeg
real::sy,sx,crow,ccol,ddx,ddy,circ,sinang,cosang
!print*,"vortex center",ccol,crow
xj=ccol-sx
yj=crow-sy
rad=sqrt(xj*xj+yj*yj)
if (rad.eq.0.) then
ddx=0.
ddy=0.
omeg=1./circ
print*,xj,yj,rad,omeg
goto 100
endif
cosang=xj/rad
sinang=yj/rad
r02=r0*r0
rad2=rad*rad
a1=r02/(2.*rad)
omeg=exp(-rad2/r02)/circ
a2=(1./circ-omeg)
a12=a1*a2
ddx=-sinang*a12
ddy=cosang*a12
if (omeg.gt.(1./circ)) print*,' bug'
100 return
end subroutine burger