反距离加群插值

懒猫的苹果

反距离加权插值的一个程序，看不太懂，哪位大神可以帮忙分析一下。。

kongfeng0824

直接贴出代码吧。这样下载估计很少有人会看

topmad

又是留作业式的提问...哎

topmad

又是留作业式的提问...哎

Aires

楼主自己先理解一下反距离加权的理论，再看程序就会轻松多了

懒猫的苹果

dimension ta(13592),txy(4,605),xn(10,10,238),ym(10,10,238)
        dimension tmn(10,10,238),mng(10,10)
        character*20 name1,name2,name3           
        open(5,file='data5')      
        write(*,*)'input name2=?'
        read(*,2555)name2
        write(*,*)'input name3=?'
        read(*,2555)name3
        write(*,*)'input name1='
        read(*,2555)name1         
        read(5,*)mn1,xk,yk,x0,x1,y0,y1,f9,rp     
        yy1=y1      
        y1=y1-y0     
        y0=0         
        n0=x0/xk+1   
        m0=y0/yk+1   
        n1=x1/xk+1   
        m1=y1/yk+1   
        n=n1-n0+1   
        m=m1-m0+1   
        mn=m*n      
        es=0.1      
        es2=es*es   
        rp2=rp*rp   
        small=1.0e-3
        small2=small*small         
        open(3,file=name3,status='new')         
        write(3,310)x0,y0,x1,y1   
        write(3,320)m,n,xk,yk      
        write(3,330)f9,rp         
310        format(1x,'x0=',f10.2,'  y0=',f10.2/1x,'x1=',f10.2,'  y1=',f10.2)
320        format(1x,'m=',i6,'  n=',i6/1x,'xk=',f6.2,'  yk=',f6.2)
330        format(1x,'f9=',f7.1,'  rp=',f6.2)   
        close(3)     
        lm=rp/yk     
        if(lm*yk.ne.rp)lm=rp/yk+1  
        if(lm.eq.0)lm=1            
        ln=rp/xk     
        if(ln*xk.ne.rp)ln=rp/xk+1  
        if(ln.eq.0)ln=1            
        mw=m/lm      
        if(mw*lm.ne.m)mw=m/lm+1   
        nw=n/ln      
        if(nw*ln.ne.n)nw=n/ln+1   
        do 10 j=1,mw+2            
        do 10 i=1,nw+2            
        write(*,*)'i=',i,'  j=',j  
10        mng(i,j)=0
        x01=x0-20   
        x11=x1+20   
        y01=y0-20   
        y11=y1+20   
        ykm=yk*float(lm)           
        xkn=xk*float(ln)           
        write(*,*)'input original name1'         
2555        format(a)              
        open(1,file=name1,status='old')
        write(0,*)'input nn=?'     
30        read(5,*)nn              
        close (5)   
        write(0,*)nn
        rewind (1)
        do 50 i=1,nn
        read(1,*)xx,yy,t
c        write(0,*)yy,xx,t         
        yy=yy1-yy   
        if(yy.le.y01.or.g.ge.y11)goto 50         
        if(xx.le.x01.or.xx.ge.x11)goto 50        
        x=xx-x0      
        y=yy-y0      
        jm=y/ykm+1   
        if(jm.lt.0.or.jm.gt.(mw+1))goto 50      
        in=x/xkn+1   
        if(in.lt.0.or.in.gt.(nw+1))goto 50      
        if(mng(in+1,jm+1).ge.200)goto 50         
        mng(in+1,jm+1)=mng(in+1,jm+1)+1         
        xn(in+1,jm+1,mng(in+1,jm+1))=x           
        ym(in+1,jm+1,mng(in+1,jm+1))=y           
        tmn(in+1,jm+1,mng(in+1,jm+1))=t         
50        continue   
60        close(1)   
        do 300 jj=2,mw+1           
        ly=(jj-2)*lm
        km=lm        
        if(jj.eq.mw+1.and.mw*lm.ne.m)km=m-mw*lm+lm            
        do 300 ii=2,nw+1           
        write(*,*)'jj  ii'         
        write(*,*)jj,ii            
        lx=(ii-2)*ln
        kn=ln        
        if(ii.eq.nw+1.and.nw*ln.ne.n)kn=n-nw*ln+ln            
        do 200 j=1,km              
        j1=ly+j      
        y=float(j1-1)*yk           
        lj=(j1-1)*n  
        do 200 i=1,kn              
        i1=lx+i      
        x=float(i1-1)*xk           
        li=lj+i1     
        ta(li)=f9   
        l=1         
        do 110 jm=jj-1,jj+1        
        do 110 in=ii-1,ii+1        
        if(mng(in,jm).eq.0)goto 110              
        do 100 k=1,mng(in,jm)      
        r2=(xn(in,jm,k)-x)**2+(ym(in,jm,k)-y)**2
        if(r2.le.rp2)then         
        if(r2.le.small2)then      
        ta(li)=tmn(in,jm,k)        
        goto 200     
        endif        
        if(l.gt.600)goto 120      
        txy(1,l)=xn(in,jm,k)      
        txy(2,l)=ym(in,jm,k)      
        txy(3,l)=tmn(in,jm,k)      
        txy(4,l)=1.0/r2            
        l=l+1        
        endif        
100        continue  
110        continue  
        l1=l-1      
120        if(l1.lt.6)goto 200     
        s1=0.0      
        s2=0.0      
        if(l1.eq.6)then            
        do 130 k=1,6
        s1=s1+txy(3,k)*txy(4,k)   
130        s2=s2+txy(4,k)         
        ta(li)=s1/s2
        goto 200     
        endif        
        do 150 k=1,6
        ic=1         
        do 140 kk=2,l1            
        if(txy(4,kk).gt.txy(4,ic))ic=kk         
140        continue  
        s1=s1+txy(3,ic)*txy(4,ic)  
        s2=s2+txy(4,ic)            
150        txy(4,ic)=1.0e-20      
        ta(li)=s1/s2
200        continue  
300        continue  
        open(2,file=name2,status='new')         
8585        format(1x,'line=',i6/(8f8.2))        
        do 5656 ii=1,m            
5656        write(2,8585)ii,(ta(ikk),ikk=(ii-1)*n+1,ii*n)      
c        write(2,8585)(ta(i),i=1,mn)            
        close(2)     
        open(9,file='wwgg',status='new')         
        do 5421 iill=1,m           
5421        write(9,6987)(ta(ikkk),ikkk=(iill-1)*n+1,iill*n)   
6987        format(1x,10f7.2)      
        stop         
        end         


fjfjkyky

{:eb513:}{:eb513:}

fjfjkyky

{:eb313:}{:eb313:}

alexclimate

谢谢分享哈。。

alexclimate

谢谢分享哈。。