diff --git a/Fortran/FUNC_median.f90 b/Fortran/FUNC_median.f90 deleted file mode 100644 index 85c7bca..0000000 --- a/Fortran/FUNC_median.f90 +++ /dev/null @@ -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 - - - \ No newline at end of file diff --git a/Fortran/FUNC_sign.f90 b/Fortran/FUNC_sign.f90 deleted file mode 100644 index 19d7669..0000000 --- a/Fortran/FUNC_sign.f90 +++ /dev/null @@ -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 - - \ No newline at end of file diff --git a/Fortran/SUB_Date_to_Num.f90 b/Fortran/SUB_Date_to_Num.f90 deleted file mode 100644 index c4f9c8f..0000000 --- a/Fortran/SUB_Date_to_Num.f90 +++ /dev/null @@ -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 \ No newline at end of file diff --git a/Fortran/SUB_FileTrans.f90 b/Fortran/SUB_FileTrans.f90 deleted file mode 100644 index aba6df6..0000000 --- a/Fortran/SUB_FileTrans.f90 +++ /dev/null @@ -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 - \ No newline at end of file diff --git a/Fortran/SUB_WAVELET.f90 b/Fortran/SUB_WAVELET.f90 deleted file mode 100644 index aa7279f..0000000 --- a/Fortran/SUB_WAVELET.f90 +++ /dev/null @@ -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 - - - - - - - - - - \ No newline at end of file diff --git a/Fortran/SUB_conv.f90 b/Fortran/SUB_conv.f90 deleted file mode 100644 index 099f3b3..0000000 --- a/Fortran/SUB_conv.f90 +++ /dev/null @@ -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 \ No newline at end of file diff --git a/Fortran/SUB_dwtLOC.f90 b/Fortran/SUB_dwtLOC.f90 deleted file mode 100644 index 1e93b36..0000000 --- a/Fortran/SUB_dwtLOC.f90 +++ /dev/null @@ -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 - - - - \ No newline at end of file diff --git a/Fortran/SUB_ebayesthresh.f90 b/Fortran/SUB_ebayesthresh.f90 deleted file mode 100644 index c718a90..0000000 --- a/Fortran/SUB_ebayesthresh.f90 +++ /dev/null @@ -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=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 deleted file mode 100644 index fcd0662..0000000 --- a/Fortran/SUB_sequence_fenduan.f90 +++ /dev/null @@ -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_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 deleted file mode 100644 index f63b62f..0000000 --- a/Fortran/SUB_upsconv.f90 +++ /dev/null @@ -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 - - \ No newline at end of file diff --git a/Fortran/SUB_weightfromdata.f90 b/Fortran/SUB_weightfromdata.f90 deleted file mode 100644 index 4e7ff26..0000000 --- a/Fortran/SUB_weightfromdata.f90 +++ /dev/null @@ -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)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