Page MenuHomeHEPForge
Files F11221259

interpol_v0-5.f90
No OneTemporary
Actions

interpol_v0-5.f90
View Options

!------------------------------------------------------------
program interpol
!------------------------------------------------------------
! This program will interpolate
! between lines of data, so that the first column is evenly spaced
!
! Minimum number of columns is 2, max is 100.
! All input must be positive.
!
! These variables need to be added by hand:
! filename: main part of filename (without file extension)
! ending: file extension i.e. will read interpol_demo.txt
! loginterpol: whether interpolation should be linear or log10
! skip: number of lines to be skipped before data starts
! ncol: number of columns of data
! power: raise all data (apart from 1st column) to this power (usually, this is 1).
!
! compile and run using:
! gfortran -C interpol_v0-5.f90 -o interpol_v0-5.exe
! ./interpol_v0-5.exe
!------------------------------------------------------------
implicit none
!------------------------------------------------------------
integer :: i,j,k,m,q,jj,y
integer :: skip,ncol
logical :: loginterpol
integer :: l_in,l_out,f_in,f_out
character(LEN=100) :: filename
double precision, allocatable :: mh_XS_in(:,:),mh_XS_out(:,:)
double precision, allocatable :: mh_log10XS_in(:,:),mh_log10XS_out(:,:)
double precision :: delta_mh,delta_mh_temp,mh_temp,mh_min,test
logical :: found_delta_mh
double precision :: delta_mh_first
character(LEN=150), allocatable :: description(:)
character(LEN=2) :: a
character(LEN=4) :: ending
integer :: power,x
!------------------------------------------------------------
!----------------------------------------------------------------
! Edit this bit!
filename='../../docs/Expt/3381_D0_h-aa_0905.3381/tableI_0905-3381_ed' !main part of filename (without file extension)
ending='.txt' !file extension i.e. will read interpol_demo.txt
loginterpol=.False. !whether interpolation should be linear or log10
skip=5 !number of lines to be skipped before data starts
ncol=3 !number of columns of data
power=1 !raise all data (apart from 1st column) to this power
!(usually, this is 1)
!E.g.
!filename='sample/interpol_demo'
!ending='.txt'
! or
!filename='../HiggsBounds_KW/Theory_tables/TeV4LHCggh-tev_ed' ; ending='.dat'
!loginterpol=.True. ; skip=0 ; ncol=2 ; power=1
! or
!filename='../HiggsBounds_KW/Expt_tables/CDFtables/CDF_H-WW_3.6fb_9023' ; ending='.txt'
!loginterpol=.False. ; skip=5 ; ncol=3 ; power=1
!----------------------------------------------------------------
if((ncol.lt.2).or.(ncol.gt.100))stop 'wrong number of columns'
if(skip.lt.0)stop 'need to set skip'
write(*,*)'starting '//trim(filename)//'...'
! open input and output files
f_in=20
f_out=f_in*10
open(f_in,file=trim(filename)//ending)
if(power.eq.1)then
open(f_out,file=trim(filename)//'_interpol'//ending)
elseif(power.eq.2)then
open(f_out,file=trim(filename)//'_sq_interpol'//ending)
else
stop'change this if you want a power that is not 1 or 2'
endif
x=getfilelength(f_in)
! find length of input file
l_in=x-skip
write(*,*)'skip=',skip
write(*,*)'l_in=',l_in
allocate(mh_XS_in(l_in,ncol))
allocate(mh_log10XS_in(l_in,ncol))
! read in input file
if(skip.ne.0)then
allocate(description(skip))
do i=1,skip
read(f_in,'(a)')description(i)
enddo
endif
do i=1,l_in
read(f_in,*)(mh_XS_in(i,j),j=1,ncol)
mh_log10XS_in(i,1)=mh_XS_in(i,1)
do j=2,ncol
mh_log10XS_in(i,j)=log10(mh_XS_in(i,j))
enddo
enddo
do i=1,l_in
do m=1,ncol
if(mh_XS_in(i,m).lt.0.0D0)stop 'negative input data'
enddo
enddo
! find smallest mass separation
delta_mh=mh_XS_in(l_in,1)
do i=1,l_in-1
if(mh_XS_in(i+1,1).le.mh_XS_in(i,1))stop 'bad input file (a)'
delta_mh_temp=mh_XS_in(i+1,1)-mh_XS_in(i,1)
if(delta_mh_temp.lt.delta_mh)delta_mh=delta_mh_temp
enddo
!need to check that delta_mh fits
!into all the other mass differences a whole number of times
y=0
found_delta_mh=.False.
delta_mh_first=delta_mh
do while((.not.found_delta_mh).and.(y.lt.5))
y=y+1
found_delta_mh=.True.
delta_mh=delta_mh_first/dble(y)
do i=1,l_in-1
delta_mh_temp=mh_XS_in(i+1,1)-mh_XS_in(i,1)
test=(delta_mh_temp/delta_mh)-dble(nint(delta_mh_temp/delta_mh))
if(abs(test).gt.1.0D-7)then
found_delta_mh=.False.
endif
enddo
enddo
if(.not.found_delta_mh)stop 'need to calculate delta_mh more carefully (find a common factor)'
write(*,*)'delta_mh=',delta_mh
! setb up output array
l_out=nint((mh_XS_in(l_in,1)-mh_XS_in(1,1))/delta_mh)+1
allocate(mh_XS_out(l_out,ncol))
allocate(mh_log10XS_out(l_out,ncol))
mh_XS_out = -1.0D0
mh_log10XS_out = -1.0D0
if(l_out.eq.l_in)then
write(*,*)' no interpolation needed for '//filename
mh_XS_out=mh_XS_in
mh_log10XS_out=mh_log10XS_in
else
! at the moment, delta_mh and mh_min are rounded to the nearest 0.1
if(abs(delta_mh-1.0D-1).gt.1.0D-5)stop'error1'
delta_mh=dble(nint(delta_mh*10.0D0)/10.0D0)
if(abs(mh_XS_in(1,1)*10.0D0-dble(nint(mh_XS_in(1,1)*10.0D0))).gt.1.0D-7)then
stop'error2'
endif
mh_min=dble(nint(mh_XS_in(1,1)*10.0D0)/10.0D0)
! fill any elements of output array which we can do without interpolation
do i=1,l_out
! fill any elements of output array which we can do without interpolation
mh_temp=dble(i-1)*delta_mh+mh_min
mh_XS_out(i,1)= mh_temp
mh_log10XS_out(i,1)= mh_temp
do j=1,l_in
if(abs(mh_XS_in(j,1)-mh_temp).lt.1.0D-7)then
do m=2,ncol
mh_XS_out(i,m)=mh_XS_in(j,m)
mh_log10XS_out(i,m)=mh_log10XS_in(j,m)
enddo
endif
enddo
enddo
! do interpolation
do i=1,l_out
if(( mh_XS_out(i,2) - (-1.0D0) ).lt.1.0D-7)then
j=i+1
do while (( mh_XS_out(j,2) - (-1.0D0) ).lt.1.0D-7)
j=j+1
enddo
do m=2,ncol
if(loginterpol)then
mh_log10XS_out(i,m)=mh_log10XS_out(i-1,m)+(mh_log10XS_out(j,m)-mh_log10XS_out(i-1,m)) &
& *(mh_XS_out(i,1)-mh_XS_out(i-1,1))/(mh_XS_out(j,1)-mh_XS_out(i-1,1))
mh_XS_out(i,m)=10.0D0**mh_log10XS_out(i,m)
else
mh_XS_out(i,m)=mh_XS_out(i-1,m)+(mh_XS_out(j,m)-mh_XS_out(i-1,m)) &
& *(mh_XS_out(i,1)-mh_XS_out(i-1,1))/(mh_XS_out(j,1)-mh_XS_out(i-1,1))
if(mh_XS_out(i,m).ge.0.0D0)mh_log10XS_out(i,m)=log10(mh_XS_out(i,m))
endif
enddo
endif
enddo
endif
if(skip.ne.0)then
do i=1,skip
! we want to leave out any comment line involving the string,
! 'interpol' since it's probably saying that the file needs interpolation,
! which is exactly what we're doing!
if(index(trim(description(i)),'interpol',.True.).eq.0)then
write(f_out,'(a)')trim(description(i))
else
write(f_out,*)
endif
enddo
deallocate(description)
endif
write(a,'(I2)')ncol-1
do i=1,l_out
do jj=2,ncol
mh_XS_out(i,jj)=mh_XS_out(i,jj)**dble(power)
enddo
write(f_out,'(1f15.2,'//trim(a)//'g15.7)')(mh_XS_out(i,j),j=1,ncol)
enddo
!********************************************
deallocate(mh_XS_in)
deallocate(mh_XS_out)
deallocate(mh_log10XS_in)
deallocate(mh_log10XS_out)
close(f_in)
close(f_out)
contains
!-------------------------------------------------------------
function getfilelength(fileid)
!-------------------------------------------------------------
implicit none
!-------------------------------------output
integer fileid
!-----------------------------------internal
integer :: getfilelength
integer :: n,ios,m
character(len=100) :: temp
!-------------------------------------------
n= 0
do !this will count the number of lines which end in a newline character
read(fileid,*,iostat=ios)
if(ios.lt.0) exit
if(ios.gt.0)then
write(*,*)'error in input file',fileid ; call flush(6)
exit
endif
n = n + 1
enddo
rewind(fileid)
m = 0 ; temp="" !this will count the number of lines in the file, including the last one,
do ! even if it doesn't end in a newline character
read(fileid,'(a)',iostat=ios)temp
if(ios.lt.0) exit
m = m + 1
enddo
rewind(fileid)
if(m.ne.n)then
write(*,*)'error in input file: file id=',fileid ; call flush(6)
stop'Error: file needs to end with a newline character (see standard output for file id)'
endif
getfilelength=n
end function getfilelength
end program interpol

File Metadata

Mime Type
text/plain
Expires
Wed, May 14, 10:09 AM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5068679
Default Alt Text
interpol_v0-5.f90 (8 KB)

Event Timeline