misaki/SUB_WAVELET
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

2

Browse Source
This commit is contained in:
misaki 2025-04-17 17:31:11 +08:00
parent fdcb8d2922
commit 4c6342ff08
12 changed files with 0 additions and 832 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
Fortran/FUNC_median.f90
View File

@ -1,33 +0,0 @@
!求数组的中位数（偶数个时取平均值）
real*8 function median(n,A)
implicit none
integer :: n
real*8 :: A(n)
real*8 :: B(n),temp
integer :: i,j
!从小到大排序
B = A
do i=1,n-1
do j=i+1,n
if(B(i)>B(j)) then
temp = B(i)
B(i) = B(j)
B(j) = temp
end if
end do
end do
!按奇偶取中位数
if(mod(n,2)==1) then
median = B(n/2+1)
else
median = (B(n/2)+B(n/2+1))/2.0
end if
end function

19
Fortran/FUNC_sign.f90
View File

@ -1,19 +0,0 @@
!符号函数
real*8 function FUNC_sign(A)
implicit none
real*8 :: A
if(isnan(A)) then
FUNC_sign = sqrt(-1.0)
elseif(A>0) then
FUNC_sign = 1.0
elseif(A<0) then
FUNC_sign = -1.0
elseif(A==0) then
FUNC_sign = 0.0
end if
end function

59
Fortran/SUB_Date_to_Num.f90
View File

@ -1,59 +0,0 @@
!日期转为数值，与matlab中Datenum函数相同
subroutine Date_to_Num(date,N,OutLine)
implicit none
integer :: date(1000,5)
integer :: N
real*8 :: OutLine(1000)
real*8 :: day1(1000),day2(1000),day3(1000),day4(1000)
real*8 :: day(1000)
integer :: i
day1 = 0
day2 = 0
day3 = 0
day4 = 0
day = 0
day1 = int((date(:,1)-1)/400) !400年个数146097
day2 = int(((date(:,1)-1)-day1*400)/100) !剩余100年的个数
day3 = int(((date(:,1)-1)-day1*400-day2*100)/4) !剩余4年的个数
day4 = int((date(:,1)-1)-day1*400-day2*100-day3*4) !剩余1年的个数
day = 366+day1*146097+day2*36524+day3*1461+day4*365
do i=1,N
select case(date(i,2))
case(2)
day(i) = day(i)+31
case(3)
day(i) = day(i)+59
case(4)
day(i) = day(i)+90
case(5)
day(i) = day(i)+120
case(6)
day(i) = day(i)+151
case(7)
day(i) = day(i)+181
case(8)
day(i) = day(i)+212
case(9)
day(i) = day(i)+243
case(10)
day(i) = day(i)+273
case(11)
day(i) = day(i)+304
case(12)
day(i) = day(i)+334
end select
if(day4(i)==3.and.(.not.(day3(i)==24.and.day2(i)/=3)).and.date(i,2)>2) then
day(i) = day(i)+1
end if
end do
day = day+date(:,3)+date(:,4)/24.0+date(:,5)/1440.0
day(N+1:) = 0.0
OutLine = day
end subroutine

78
Fortran/SUB_FileTrans.f90
View File

@ -1,78 +0,0 @@
!按格式读取文件，将csv文件按字符串形式分列读取，对时间格式进行处理
!csv文件首列为时间，后续为数值，参考智慧大坝环境文件
!时间格式支持年月日时分和年月日，年月日以/隔开，时分以:隔开
!若数值有多列，字符串长度适度增加
subroutine FileTrans(M,FileLen,FileName,N,OutLine)
implicit none
common /tips/err
integer :: err
integer :: M !文件列数
integer :: FileLen
character(len=FileLen) :: FileName
real*8 :: OutLine(1000,M) !时间数值（精确到小时）超过8位
character(len=200) :: cha
character(len=200) :: cdate,cnum
integer :: date(1000,5)
real*8 :: num(1000,M-1)
integer :: i,stat,j
integer :: datenum !日期格式数量
integer :: N !数列长度
N = 0
datenum = 0
date = 0
num = 0
OutLine = 0
open(100,file=FileName,status='old')
read(100,'(a)')cha !第一行为表头
do i=1,40000
read(100,'(a)',iostat=stat)cha
if(cha(len_trim(cha):len_trim(cha))==',') then
err = 1 !数据缺失(末行)
goto 100
end if
if(stat/=0) exit
N = N+1
if(index(cha(1:len_trim(cha)),':')>0) then
datenum = 5
else
datenum = 3
end if
if(index(cha,',,')>0.or.index(cha,'/')==0) then
err = 1 !数据缺失（首行或中间行）
goto 100
end if
do j=1,len_trim(cha)
if(index('-/:',cha(j:j))>0) cha(j:j)=' '
if(cha(j:j)==',') then
cdate = cha(1:j-1)
cnum = cha(j+1:)
exit
end if
end do
read(cdate,*)date(i,1:datenum)
read(cnum,*)num(i,:)
end do
close(100)
!当计算月份时，需要获取初始时间，存放在序列最后一行，所以此时文件数据不能超过999行
date(N+1,1) = date(1,1)
date(N+1,2:3) = 1
OutLine(:,2:M) = num(:,1:M-1)
call Date_to_num(date,N,OutLine(:,1))
100 end subroutine

82
Fortran/SUB_WAVELET.f90
View File

@ -1,82 +0,0 @@
SUBROUTINE SUB_WAVELET(n,&
x,&
xden)&
bind(C, name="SUB_WAVELET")
!DEC$ ATTRIBUTES DLLEXPORT::SUB_WAVELET
implicit none
integer :: n
real*8 :: x(n),xden(n)
integer :: level,sz,temp,num,temp_n,p,s,i
real*8 :: Lo_D(8),Hi_D(8),Lo_R(8),Hi_R(8),a(1000),d(1000),temp_xfinal(1000),normfac
real*8,allocatable :: origcfs(:,:),cd(:,:)
integer,allocatable :: temp_sx(:)
real*8 :: temp_x(1000),temp_x2(1000),median,vscale
xden = 0
level = int(log(n/7.0)/log(2.0))
Lo_D = [-0.0757657147892733,-0.0296355276459985,0.497618667632015,0.803738751805916,0.297857795605277,-0.0992195435768472,-0.0126039672620378,0.0322231006040427]
Hi_D = [-0.0322231006040427,-0.0126039672620378,0.0992195435768472,0.297857795605277,-0.803738751805916,0.497618667632015,0.0296355276459985,-0.0757657147892733]
Lo_R = [0.0322231006040427,-0.0126039672620378,-0.0992195435768472,0.297857795605277,0.803738751805916,0.497618667632015,-0.0296355276459985,-0.0757657147892733]
Hi_R = [-0.0757657147892733,0.0296355276459985,0.497618667632015,-0.803738751805916,0.297857795605277,0.0992195435768472,-0.0126039672620378,-0.0322231006040427]
call dwtLOC(n,x,Lo_D,Hi_D,0,2,a,d,sz)
allocate(origcfs(1000,level+1),temp_sx(level+2))
temp_xfinal = 0
origcfs = 0
temp_sx = 0
temp_sx(1) = n !şón-1¸öĘýĘÇorigcfsľÄĘýžÝ¸öĘý
do i=level+1,2,-1
if(i==level+1) then
call dwtLOC(n,x,Lo_D,Hi_D,0,2,a,d,sz)
temp_xfinal = 0
temp_xfinal(1:sz) = a(1:sz)
else
call dwtLOC(sz,temp_xfinal,Lo_D,Hi_D,0,2,a,d,temp)
temp_xfinal = 0
sz = temp
temp_xfinal(1:sz) = a(1:sz)
end if
origcfs(1:sz,level+2-i) = d(1:sz)
temp_sx(level+2-i+1) = sz
if(i==2) temp_sx(level+2) = sz
end do
origcfs(1:sz,level+1) = temp_xfinal(1:sz)
normfac = 1.48260221850560
vscale = normfac*median(temp_sx(2),abs(origcfs(1:temp_sx(2),1)))
allocate(cd(1000,level))
cd = 0
do i=1,level
call ebayesthresh(temp_sx(i+1),origcfs(1:temp_sx(i+1),i),vscale,cd(1:temp_sx(i+1),i))
end do
temp_x = origcfs(:,level+1)
temp_n = temp_sx(level+2)
do i=level,1,-1
p = level+2-i
s = temp_sx(level+2-p)
call upsconv(temp_n,temp_x,Lo_R,s,temp_x2)
call upsconv(temp_sx(i+1),cd(1:temp_sx(i+1),i),Hi_R,s,temp_x)
temp_x = temp_x+temp_x2
temp_n = s
end do
xden = temp_x(1:n)
end subroutine

18
Fortran/SUB_conv.f90
View File

@ -1,18 +0,0 @@
!蕨좆얩생
SUBROUTINE conv(m,x,n,y,z)
implicit none
integer :: m,n
real*8 :: x(m),y(n),z(m+n-1)
integer :: i,j
z = 0
do i=1,m
do j=1,n
z(i+j-1) = z(i+j-1)+x(i)*y(j)
end do
end do
end subroutine

65
Fortran/SUB_dwtLOC.f90
View File

@ -1,65 +0,0 @@
SUBROUTINE dwtLOC(n,x,LoD,HiD,perFLAG,first,a,d,sz)
implicit none
integer :: n,perFLAG,first,sz
real*8 :: x(n),LoD(8),HiD(8),a(1000),d(1000)
integer :: dCol,lenEXT,lenKEPT,i,j
integer,allocatable :: idxCOL(:)
real*8,allocatable :: y(:),c1(:),c2(:)
a = 0
d = 0
sz = 0
dCol = 7
if(perFLAG==0) then
lenEXT = 7
lenKEPT = n+7
else
lenEXT = 4
lenKEPT = 2*int(ceiling(n/2.0))
end if
sz = int(lenKEPT/2.0)
allocate(idxCOL(sz))
do i=1,sz
idxCOL(i) = first+dCol+2*(i-1)
end do
allocate(y(n+2*lenEXT))
y(1:lenEXT) = x(lenEXT:-1:1)
do i=1,lenEXT
y(i) = x(lenEXT+1-i)
end do
do i=lenEXT+1,lenEXT+n
y(i) = x(i-lenEXT)
end do
do i=lenEXT+n+1,n+2*lenEXT
y(i) = x(2*n+lenEXT+1-i)
end do
!y(lenEXT+1:lenEXT+n) = x
!y(lenEXT+n+1:n+2*lenEXT) = x(n:-1:n-lenEXT+1)
allocate(c1(n+2*lenEXT+7),c2(n+2*lenEXT+7))
c1 = 0
c2 = 0
do i=1,n+2*lenEXT
do j = 1,8
c1(i+j-1) = c1(i+j-1)+y(i)*LoD(j)
c2(i+j-1) = c2(i+j-1)+y(i)*HiD(j)
end do
end do
do i=1,sz
a(i) = c1(idxCOL(i))
d(i) = c2(idxCOL(i))
end do
end subroutine

38
Fortran/SUB_ebayesthresh.f90
View File

@ -1,38 +0,0 @@
SUBROUTINE ebayesthresh(m,x,stdest,muhat)
implicit none
integer :: m
real*8 :: x(m),stdest,muhat(m)
real*8 :: minstd,temp_stdest(m),tempx(m),weight,pi
integer :: maxiter
maxiter = 50
minstd = 1.0e-9
pi = 3.1415926
if(stdest<minstd) then
temp_stdest = minstd
else
temp_stdest = stdest
end if
tempx = x/temp_stdest
call weightfromdata(m,tempx,30,weight)
call postmedcauchy(m,tempx,weight,maxiter,muhat)
muhat = muhat*temp_stdest
end subroutine

86
Fortran/SUB_postmedcauchy.f90
View File

@ -1,86 +0,0 @@
SUBROUTINE postmedcauchy(m,x,weight,maxiter,muhat)
implicit none
integer :: m,maxiter
real*8 :: x(m),weight,muhat(m)
real*8 :: magdata(m),magdatatmp(m),weight2(m),hi(m),Tol,conTol,lo(m),midpoint(m),y(m),pi,Xvec(m),normconstants,fx(m),delta(maxiter)
real*8 :: z(m),y2(m),y3(m),fx2(m),fmidpoint(m),zeromd(m),mu,sigma,yr(m),FUNC_sign
integer :: i,idx(m),numiter
muhat = 0
magdata = abs(x)
magdatatmp = magdata
idx = 1
delta = 0
do i=1,m
if(magdata(i)>=20) then
idx(i)=0
magdata(i) = sqrt(-1.0)
end if
end do
weight2 = weight
lo = 0
hi = maxval(magdata)
Tol = 1e-9
numiter = 0
conTol = 1e30
mu = 0
sigma = 1
pi = 3.1415926
zeromd = 0
do while(conTol>Tol)
numiter = numiter +1
midpoint = (lo+hi)/2.0
y = magdata-midpoint
Xvec = abs((y-mu)**2)
normconstants = 1.0/(sqrt(2.0*pi)*sigma)
fx = normconstants*exp(-Xvec/(2.0*sigma))
z = (y-mu)/sigma
y2 = 1.0/2.0*erfc(-z/sqrt(2.0))
z = (-midpoint-mu)/sigma
y3 = 1.0/2.0*erfc(-z/sqrt(2.0))
Xvec = abs((midpoint-mu)**2)
fx2 = normconstants*exp(-Xvec/(2.0*sigma))
yr = y2 - magdata*fx + ((magdata*midpoint-1)*fx*y3/fx2)
y3 = magdata**2.0*(1.0/weight2-1.0)-1.0
y2 = 1+exp(-magdata**2/2.0)*y3
fmidpoint = y2/2.0-yr
do i=1,m
if(fmidpoint(i) <= zeromd(i)) then
lo(i) = midpoint(i)
else
hi(i) = midpoint(i)
end if
end do
delta(numiter) = maxval(abs(hi-lo))
conTol = delta(numiter)
if(numiter > maxiter) exit
end do
muhat = (lo+hi)/2.0
do i=1,m
if(idx(i)==0) then
muhat(i) = magdatatmp(i)-2.0/magdatatmp(i)
end if
if(muhat(i)<1e-7) muhat(i)=0
muhat(i) = FUNC_sign(x(i))*muhat(i)
if(abs(muhat(i)) > abs(x(i))) muhat(i) = x(i)
end do
end subroutine

239
Fortran/SUB_sequence_fenduan.f90
View File

@ -1,239 +0,0 @@
SUBROUTINE sequence_fenduan(n,date,deformation,judge,num,trendvalue,dn)
implicit none
integer :: n,judge,num(1000),dn
real*8 :: date(n),deformation(n),trendvalue(2000)
real*8 :: X_Time(n),xden(n)
real*8,allocatable :: B(:,:),X(:,:),X1(:),X2(:),data_set(:),data_setcheck(:),time_set(:),tempdata(:),Y(:,:),Y_check(:,:)
real*8,allocatable :: TrendValue2(:,:),error_data(:),error_data2(:),error_datackeck(:),FF(:),AA(:)
real*8 :: intercept
integer :: xn,i,VV,VV1,VV2,tempn,location,loca,loca2,loca_check
integer,allocatable :: inmodel(:,:),datanum(:),tempnum(:)
real*8 :: FF_check,finv(1000,3),X_check
integer :: sn !总数据量
num = 0
trendvalue = 0
dn = 0
X_Time = (date-(date(1)-1))/100.0
!小波函数去噪
call wdenoise(n,deformation,xden)
sn = 3*(n-2)
dn = n-2
allocate(data_set(sn),data_setcheck(sn),time_set(sn),datanum(dn))
datanum = 3
do i=1,n-2
data_set((i-1)*3+1:i*3) = Deformation(i:i+2)
data_setcheck((i-1)*3+1:i*3) = xden(i:i+2)
time_set((i-1)*3+1:i*3) = X_Time(i:i+2)
end do
!读取逆累积分布函数数据
open(20,file='finv.txt',status='old')
do i=1,1000
read(20,*)finv(i,:)
end do
close(20)
!逐步回归拟合
do while(.true.)
allocate(error_data(dn-1),error_data2(dn-1),error_datackeck(dn-1),AA(dn-1))
do i=1,dn-1
tempn = datanum(i)+datanum(i+1)-3+1
allocate(X1(tempn),X2(tempn),Y(tempn,1),Y_check(tempn,1))
X1(1:datanum(i)) = time_set(sum(datanum(1:i-1))+1:sum(datanum(1:i)))
X1(datanum(i)+1:tempn) = time_set(sum(datanum(1:i))+3:sum(datanum(1:i+1)))
X2 = X1**2
if(judge==1) then
allocate(X(tempn,1))
X(:,1) = X1
xn = 1
elseif(judge==2) then
allocate(X(tempn,2))
X(:,1) = X1
X(:,2) = X2
xn = 2
elseif(judge==0) then
allocate(X(tempn,1))
X = 1
xn = 1
end if
Y(1:datanum(i),1) = data_set(sum(datanum(1:i-1))+1:sum(datanum(1:i)))
Y(datanum(i)+1:tempn,1) = data_set(sum(datanum(1:i))+3:sum(datanum(1:i+1)))
Y_check(1:datanum(i),1) = data_setcheck(sum(datanum(1:i-1))+1:sum(datanum(1:i)))
Y_check(datanum(i)+1:tempn,1) = data_setcheck(sum(datanum(1:i))+3:sum(datanum(1:i+1)))
allocate(B(1,xn),inmodel(1,xn))
B = 0
inmodel = 0
call STEPWISEFIT(tempn,xn,B,inmodel,intercept,X,Y,0.05d0,0.1d0)
!sb(i) = intercept
!if(judge==2) then
! sb(2:3,i) = transpose(B*inmodel)
!else
! sb(2:3,i) = B(1,1)*inmodel(1,1)
!end if
allocate(TrendValue2(tempn,1))
TrendValue2 = intercept+matmul(X,transpose(B)*inmodel)
error_data(i) = sum((TrendValue2(:,1)-Y(:,1))**2)
error_data2(i) = (1.0/(4-1-judge))*error_data(i)
error_datackeck(i) = (1.0/(4-1))*sum((Y(:,1)-Y_check(:,1))**2)
deallocate(X,B,inmodel,X1,X2,Y,TrendValue2,Y_check)
end do
!自适应参数结束循环
allocate(FF(dn-1))
FF = error_data2/error_datackeck
FF_check = minval(FF)
loca_check = minloc(FF,1)
deallocate(FF)
VV = datanum(loca_check)+datanum(loca_check+1)-2
VV1 = VV-1-judge
VV2 = VV-1
!逆累计分布函数，换成分类赋值
if(judge==0) then
if(VV<=1000) then
X_check = finv(VV,1)
else
X_check = finv(1000,1)+(VV-1000)*(finv(1000,1)-finv(999,1))
end if
elseif(judge==1) then
if(VV<=1000) then
X_check = finv(VV,2)
else
X_check = finv(1000,2)+(VV-1000)*(finv(1000,2)-finv(999,2))
end if
elseif(judge==2) then
if(VV<=1000) then
X_check = finv(VV,3)
else
X_check = finv(1000,3)+(VV-1000)*(finv(1000,3)-finv(999,3))
end if
end if
!判断残差平方和是否小于给定阈值，如果小于，合并两相邻段，否则退出循环
if(FF_check<x_check) then
call unique(dn-1,error_data,tempn,AA)
do i=1,dn-1
if(error_data(i)==AA(1)) then
loca = i
exit
end if
end do
if(tempn>1) then
do i=1,dn-1
if(error_data(i)==AA(2)) then
loca2 = i
exit
end if
end do
else
exit
end if
!误差最小时为最后一组数据时，取误差第二小的数据与后面一组数据合并
if(loca==dn) then
location = loca2
else
location = loca
end if
!数据重组
allocate(tempdata(sn-2),tempnum(dn-1))
tempn = sum(datanum(1:location-1))+datanum(location)-2
tempdata(1:tempn) = data_set(1:tempn)
tempdata(tempn+1:sn-2) = data_set(sum(datanum(1:location))+1:sn)
deallocate(data_set)
allocate(data_set(sn-2))
data_set = tempdata
tempdata(1:tempn) = data_setcheck(1:tempn)
tempdata(tempn+1:sn-2) = data_setcheck(sum(datanum(1:location))+1:sn)
deallocate(data_setcheck)
allocate(data_setcheck(sn-2))
data_setcheck = tempdata
tempdata(1:tempn) = time_set(1:tempn)
tempdata(tempn+1:sn-2) = time_set(sum(datanum(1:location))+1:sn)
deallocate(time_set)
allocate(time_set(sn-2))
time_set = tempdata
tempnum(1:location-1) = datanum(1:location-1)
tempnum(location) = datanum(location)-2+datanum(location+1)
tempnum(location+1:dn-1) = datanum(location+2:dn)
deallocate(datanum)
allocate(datanum(dn-1))
datanum = tempnum
deallocate(tempdata,tempnum,error_data,error_data2,error_datackeck,AA)
sn = sn-2
dn = dn-1
else
deallocate(error_data,error_data2,error_datackeck,AA)
exit
end if
end do
do i=1,dn
allocate(X1(datanum(i)),X2(datanum(i)),Y(datanum(i),1))
X1 = time_set(sum(datanum(1:i-1))+1:sum(datanum(1:i)))
X2 = X1**2
if(judge==1) then
allocate(X(datanum(i),1))
X(:,1) = X1
xn = 1
elseif(judge==2) then
allocate(X(datanum(i),2))
X(:,1) = X1
X(:,2) = X2
xn = 2
elseif(judge==0) then
allocate(X(datanum(i),1))
X = 1
xn = 1
end if
Y(:,1) = data_set(sum(datanum(1:i-1))+1:sum(datanum(1:i)))
allocate(B(1,xn),inmodel(1,xn))
B = 0
inmodel = 0
call STEPWISEFIT(datanum(i),xn,B,inmodel,intercept,X,Y,0.05d0,0.1d0)
allocate(TrendValue2(datanum(i),1))
TrendValue2 = intercept+matmul(X,transpose(B)*inmodel)
TrendValue(sum(datanum(1:i-1))+1:sum(datanum(1:i))) = TrendValue2(:,1)
deallocate(TrendValue2,X1,X2,X,Y,B,inmodel)
end do
num(1:dn) = datanum(1:dn)
deallocate(data_set,data_setcheck,time_set,datanum)
end subroutine

30
Fortran/SUB_upsconv.f90
View File

@ -1,30 +0,0 @@
SUBROUTINE upsconv(n,x,f,s,y)
implicit none
integer :: n,s
real*8 :: x(n),f(8),y(1000)
integer :: sx2,first,last,lenKept,i
real*8 :: tempy(1000),d
y = 0
sx2 = 2*n
first = 0
last = 0
lenKept = s
do i=1,sx2-1,2
y(i) = x((i+1)/2)
end do
call conv(sx2-1,y(1:sx2-1),8,f,tempy(1:sx2+6))
if(lenKept>sx2+6) lenKept = sx2+6
d = (sx2+6-lenKept)/2.0
first = 1+int(floor(d))
last = sx2+6-int(ceiling(d))
y(1:s) = tempy(first:last)
end subroutine

85
Fortran/SUB_weightfromdata.f90
View File

@ -1,85 +0,0 @@
SUBROUTINE weightfromdata(m,x,maxiter,weight)
implicit none
integer :: m,maxiter
real*8 :: x(m),weight
real*8 :: thr,mu,sigma,Xvec,normconstants,fx,z,bfx,wlo
real*8 :: pi,phix(m),dnorm0,shi,slo,temp(m),beta(m),initialwlo,tmpweight,temp_deltaweight(100),wmid,smid,whi,conTol,wtol,stol
integer :: ii,tnum,i
pi = 3.1415926
thr = sqrt(2.0*log(real(m)))
mu = 0
sigma = 1
Xvec = abs((thr-mu)**2)
normconstants = 1.0/(sqrt(2.0*pi)*sigma)
fx = normconstants*exp(-Xvec/(2.0*sigma))
z = (thr-mu)/sigma
bfx = 1.0/2.0*erfc(-z/sqrt(2.0))
wlo = 1 + (bfx - thr*fx - 1.0/2)/(sqrt(pi/2.0)*fx*thr**2)
if(isnan(wlo)) wlo = 1
wlo = 1.0/wlo
phix = abs((x-mu)**2)
phix = normconstants*exp(-phix/(2.0*sigma))
dnorm0 = 0.398942280401433
do i=1,m
if(abs(phix(i) - dnorm0)>5.55111512312578e-14) then
beta(i) = (dnorm0/phix(i)-1.0)/x(i)**2-1.0
else
beta(i) = -0.5
end if
end do
shi = sum(beta/(1.0 + beta))
if(shi>=0) then
weight = 1
goto 10
end if
temp = wlo*beta
slo = sum(beta/(1.0+temp))
if(slo<=0) then
tmpweight = wlo
initialwlo = wlo
end if
conTol = 1.0e30
wtol = 2.22044604925031e-14
stol = 1e-7
ii = 1
temp_deltaweight = 0
temp_deltaweight(1) = 1.0-wlo
tnum = 1
whi = 1
do while(conTol>wtol)
wmid = sqrt(whi*wlo)
smid = sum(beta/(1.0+wmid*beta))
if(abs(smid)<stol) then
tmpweight = wmid
weight = tmpweight
goto 10
end if
if(smid>0) wlo = wmid
if(smid<0) whi = wmid
tnum = tnum+1
temp_deltaweight(tnum) = whi-wlo
conTol = abs(temp_deltaweight(ii+1) - temp_deltaweight(ii))
ii = ii+1
if(ii>maxiter) exit
end do
tmpweight = sqrt(wlo*whi)
if(slo <= 0) tmpweight = initialwlo
weight = tmpweight
10 end subroutine