24小时热门版块排行榜    

查看: 723  |  回复: 1

dailjiao0706

新虫 (初入文坛)

[求助] fortran中如何设计一个程序对内插法函数进行调用·不会编外部的代码,求高数帮忙已有1人参与

x轴是EK数值(放在文件中,有100项),y轴是fh数值(放在文件中,有100项),因为x,y间隔不均匀所以采用内插法函数将间隔变均匀,其中内插法函数已有,问题是编一个程序,将这个内插法函数成功调用最终得到均匀的xy数值。已知EK值及fh值单行放于txt文件中。下面这是内插法函数。
IMPLICIT NONE
real(kind=8), PARAMETER :: pi = 3.1415926535897932    !pi
real(kind=8), PARAMETER :: c = 2.99792458e08          !Speed of light (m/s)
real(kind=8), PARAMETER :: e0=1.0e07/(4*pi*c**2)      !vacuum permittivity
real(kind=8), PARAMETER :: er=1                                          !relative permittivity
real(kind=8), PARAMETER :: u0=pi*4e-07                              !vacuum permeability
real(kind=8), PARAMETER :: mp=1.6726231e-27           !Proton rest mass(kg)
real(kind=8), PARAMETER :: me=9.1093897e-31           !Electron rest mass(kg)
real(kind=8), PARAMETER :: e=1.60217733e-19           !Elementary charge (C)
END MODULE  my_Constants
!====================================================================================
!
REAL(kind=8) FUNCTION interpolation_c(x,xa,ya,NI) RESULT(retval)
implicit none
real(kind=8),intent(IN):: x
integer,intent(IN)::NI
real(kind=8),intent(IN),DIMENSION (NI):: xa,ya
!real,intent(OUT):: retval
integer:: cur,top,bot
real(kind=8):: numer,denom
logical::dec=0
!
        cur=floor(real(NI)/2)
        top=NI
        bot=1
!
top_if:  if (top>1) then
        dec=(xa(1)>xa(top))
        not_if:if(((.NOT. dec) .AND. (x<xa(1))) .OR.(dec .and. (x>xa(1)))) then
                top=1
        else
                do while(top>bot)
                        log1_if: if(((.NOT. dec) .AND. (x>xa(cur))) .OR.(dec .and. (x<xa(cur)))) then
                                bot=cur
                                cur=floor(real(cur+top)/2)
                        else
                                top=cur
                                cur=floor(real(bot+cur)/2)
                        end if log1_if
                        log2_if: if(top>(bot+1)) then
                                cur=floor(real(bot+top)/2)
                        else
                                bot=top-1
                                exit
                        end if log2_if
                end do
        end if not_if
end if top_if
!----------------------------------
n_if: if (NI <= 1) then
        retval= ya(bot)
end if n_if
eq_if: if (xa(top) == xa(bot)) then
        eq1_if: if (bot == 0) then
                retval= ya(bot)
        end if eq1_if
        retval=ya(top)
end if eq_if
        numer = ( (ya(top) - ya(bot)) * (x - xa(bot)) )
        denom = (xa(top) - xa(bot))
        retval = ya(bot) + ( numer / denom )
END FUNCTION interpolation_c
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fish.yfyh

铜虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +1
dailjiao0706: 金币+5, ★★★很有帮助 2013-12-17 16:53:23
在写程序接口的时候,你需要确认input是什么?output是什么?然后就好办了。
比如你的input是文件名(当然对应文件已经包含数据),而输出就是你想要的数据xy。
那么接口子程序可以写为
subroutine my_interpolation_c(filename, xy)
然后在这个子程序里面,再去调用 interpolation_c这个函数即可.
2楼2013-12-17 11:28:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 dailjiao0706 的主题更新
信息提示
请填处理意见