diff --git a/Fortran/FUNC_median.f90 b/Fortran/FUNC_median.f90 new file mode 100644 index 0000000..85c7bca --- /dev/null +++ b/Fortran/FUNC_median.f90 @@ -0,0 +1,33 @@ + + !求数组的中位数（偶数个时取平均值） + 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 + + + \ No newline at end of file diff --git a/Fortran/FUNC_sign.f90 b/Fortran/FUNC_sign.f90 new file mode 100644 index 0000000..19d7669 --- /dev/null +++ b/Fortran/FUNC_sign.f90 @@ -0,0 +1,19 @@ + + !符号函数 + 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 + + \ No newline at end of file diff --git a/Fortran/SUB_Date_to_Num.f90 b/Fortran/SUB_Date_to_Num.f90 new file mode 100644 index 0000000..c4f9c8f --- /dev/null +++ b/Fortran/SUB_Date_to_Num.f90 @@ -0,0 +1,59 @@ + + !日期转为数值，与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 \ No newline at end of file diff --git a/Fortran/SUB_FileTrans.f90 b/Fortran/SUB_FileTrans.f90 new file mode 100644 index 0000000..aba6df6 --- /dev/null +++ b/Fortran/SUB_FileTrans.f90 @@ -0,0 +1,78 @@ + + !按格式读取文件，将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 + \ No newline at end of file diff --git a/Fortran/SUB_WAVELET.f90 b/Fortran/SUB_WAVELET.f90 new file mode 100644 index 0000000..aa7279f --- /dev/null +++ b/Fortran/SUB_WAVELET.f90 @@ -0,0 +1,82 @@ + + 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 + + + + + + + + + + \ No newline at end of file diff --git a/Fortran/SUB_conv.f90 b/Fortran/SUB_conv.f90 new file mode 100644 index 0000000..099f3b3 --- /dev/null +++ b/Fortran/SUB_conv.f90 @@ -0,0 +1,18 @@ + + !向量卷积 + 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 \ No newline at end of file diff --git a/Fortran/SUB_dwtLOC.f90 b/Fortran/SUB_dwtLOC.f90 new file mode 100644 index 0000000..1e93b36 --- /dev/null +++ b/Fortran/SUB_dwtLOC.f90 @@ -0,0 +1,65 @@ + + 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 + + + + \ No newline at end of file diff --git a/Fortran/SUB_ebayesthresh.f90 b/Fortran/SUB_ebayesthresh.f90 new file mode 100644 index 0000000..c718a90 --- /dev/null +++ b/Fortran/SUB_ebayesthresh.f90 @@ -0,0 +1,38 @@ + + 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=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 + \ No newline at end of file diff --git a/Fortran/SUB_sequence_fenduan.f90 b/Fortran/SUB_sequence_fenduan.f90 new file mode 100644 index 0000000..fcd0662 --- /dev/null +++ b/Fortran/SUB_sequence_fenduan.f90 @@ -0,0 +1,239 @@ + + 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_check1) 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 + + + + + + + + + + + + + + \ No newline at end of file diff --git a/Fortran/SUB_upsconv.f90 b/Fortran/SUB_upsconv.f90 new file mode 100644 index 0000000..f63b62f --- /dev/null +++ b/Fortran/SUB_upsconv.f90 @@ -0,0 +1,30 @@ + + 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 + + \ No newline at end of file diff --git a/Fortran/SUB_weightfromdata.f90 b/Fortran/SUB_weightfromdata.f90 new file mode 100644 index 0000000..4e7ff26 --- /dev/null +++ b/Fortran/SUB_weightfromdata.f90 @@ -0,0 +1,85 @@ + + 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)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 \ No newline at end of file