blog
美国广域增强系统的电离层算法
WAAS（Wide Area Augmentation System）是由美国联邦航空局开发建立的用于空中导航的一个系统， 该系统主要是通过解决广域差分GPS的数据通信问题从而提高全球定位系统的精度和可用性。
美国联邦航空主管机关（FAA）及交通部是WAAS程序的发展者, 并应用于精准的航空器入场之用。 一般的GPS并不具有接收FAA的系统能力，WAAS是分析电离层的干扰所导致的讯号时差，及计算卫星轨道误差值来校正GPS讯号的， 并提供每一个卫星状况的重要数据，虽然WAAS还没有被航空界认证，这个系统还是可供给一般民众使用，如船舶及平常的GPS用户。
WAAS（图片来自网络）
WAAS由约25个地面参考站台所组成，位置横跨整个美国。它们监看着卫星的数据。 其中两个主要的站台，位于两岸，从各参考站台收集校正讯号并建立校正信息。 这些信息根据卫星的轨道误差，卫星上的电子钟误差，以及由于大气及电离层所造成的讯号延迟。 这些数据经计算校正之后再经由一个或两个地球同步卫星传播，或是轨道固定在赤道上空的卫星。 这些信息的结构是完全兼容于基本的GPS讯号，这意味着任何其有WAAS启始的GPS接收机将能解读这种讯号。
其中电离层改正是WAAS系统的一项重要内容，这里编写了WAAS系统电离层算法，仅供参考。
定义模块
module iono_grid
!----------------------------------------module coment
! Version : V1.0
! Coded by : song.yz
! Date : 2014-03-04 16:20:31 *星期二*
!-----------------------------------------------------
! Description :
! data structureof Grid point
!-----------------------------------------------------
implicit real*8(a-h,o-z)
save
!-----------data structure of Grid piont
type gridpoint !structure of one point
real*8::xlong !geographic longitude (unit: degree)
real*8::xlat !geographic latitude (unit: degree)
end type gridpoint
type (gridpoint):: grid(202,0:10)
!------------------------------------
!-------------------------------------interference for outside
type bandgrid
type (gridpoint):: apoint(202)
end type bandgrid
type bandsgrid
type (bandgrid) ::bandpoint(0:10)
end type bandsgrid
!-------------------------------------
end module iono_grid
subroutine set_grid (grid_out)
!---------------------------------subroutine comment
! Purpose : set the Grid structure matrix
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-04 16:26:47
!
!-----------------------------------------------------
! Input Parameters :
!
!
! Output Parameters :
!
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use iono_grid
implicit real*8(a-h,o-z)
type (bandsgrid) :: grid_out
character*20::ckey
open(unit=1001,file='./data/grid.dat',status='old',iostat=ierror)
!write(*,*) 'open grid.data correct ?' ,ierror
do
read(1001,*,iostat=ioend) ckey
if (ioend /= 0) exit
! exit the loop if reach the end of the file
if (ckey(1:5)=='begin') exit
! exit if get the key word of 'begin'
end do
! get the longitude and latitude by reading the file
do
read(1001,*,iostat=ioend) iband,iindex,xlong,xlat
if (ioend /= 0) exit
!exit if reach to the end of the file
!grid(0:10,202)
grid(iindex,iband)%xlong=xlong
grid(iindex,iband)%xlat=xlat
! write(*,*) grid(iindex,iband)%xlong
! write(*,*) grid(iindex,iband)%xlat
!
end do
!rewind(1001)
close(1001)
!-------------------------output the Grid info.
do i=0,10
do j=1,202
grid_out%bandpoint(i)%apoint(j)%xlong=grid(j,i)%xlong
grid_out%bandpoint(i)%apoint(j)%xlat =grid(j,i)%xlat
end do
end do
!-----------------------------------------------
end subroutine set_grid
定义结构体
module def_fi
!----------------------------------------module coment
! Version : V1.0
! Coded by : song.yz
! Date : 2014-04-01 09:40:40 *星期二*
!-----------------------------------------------------
! Description :
! data structure of iono_grid
! Post Script :
! 1.
! 2.
!
!-----------------------------------------------------
! Contains :
! *.
! *.
!-----------------------------------------------------
implicit real*8(a-h,o-z)
save
integer,parameter:: MAX_PIERCE_NO=50000
type ipptype
sequence
integer*1:: iflagAB
integer*1 :: ierr
integer*1 :: SysType
integer*2 :: week
integer*4 :: sow
real*8 :: phi
real*8 :: lam
real*8 :: IncVtec
real*8 :: RcIono
real*8 :: Elevation
real*8 :: Azimuth
integer*1 :: stationNo
integer*1 :: rcvNo
integer*1 :: satelliteNo
integer*1 :: choseWN
integer*1 :: choseB
real*8::cVtec
end type ipptype
type ipplist
sequence
integer*1 :: FlagAB
integer*4 :: n_pierce !实际穿刺点个数
type(ipptype)::pierce(MAX_PIERCE_NO)
end type ipplist
type (ipplist) :: ipp
!所有穿刺点数据结构
!input ---------
! out put-----------
!----------
! the following Data struct are defined by song.yezhi
type GridData
integer*4::GridID
real*8::Delay
integer*4::GiveI
! real*8::Give !defined by song.yz
end type GridData
type BandData
integer*1::BandID
integer*1::Mask(210) !一个带中的掩码，当Mask(i)=1表示算出了格网点延迟，
!等于0表示未算出格网点延迟
type(GridData)::Grid(210)
end type BandData
type IGPlist
integer*2::Week
integer*4::Sow
integer*1::BandNum
integer*1::SysType
type(BandData)::GridBand(0:10)
end type IGPlist
type (IGPlist) :: iGP
! 所有格网点数据结构
end module def_fi
subroutine ionogrid(iGP_p,ipp_p,week,sow,Grid_p)
!---------------------------------subroutine comment
! Purpose : 由穿刺点计算格网点延迟及GIVEI
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-12 15:26:21
! 周和秒由外部传入
!
! *** 说明：需要配置文件 gird.dat 不能删除
! 该函数为计算主函数
!-----------------------------------------------------
! Input Parameters :
! ipp-----------------穿刺点结构体
! week----------------格网点周
! sow-----------------格网点秒
! week 和 sow 由外部传入
! Grid_p-----------------格网结构
! Output Parameters :
! IGP-----------------格网点结构体
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use iono_grid
use def_fi
implicit real*8(a-h,o-z)
real*8,parameter:: Factor = 0.1651d0
!基于B1频点 TECU 到 长度单位（米） 转换系数
integer *2 week
integer *4 sow
type (Ipplist) :: iPP_p
type (IGPlist) :: iGP_p
! type (gridpoint) :: Grid_p
! type (bandsgrid) :: grid_p(0:10)
type (bandsgrid) :: grid_p
!type (gridpoint):: grid_song(202,0:10)
! Grid=Grid_p
!该部分由FORTRAN读入，并放在公共内存，这里传回给Grid
!----------------------------------------------------
do i=0,10
do j=1,202
Grid(j,i)%xlong=Grid_p%bandpoint(i)%apoint(j)%xlong
Grid(j,i)%xlat =Grid_p%bandpoint(i)%apoint(j)%xlat
end do
end do
!-----------------------------------------------------
!------------------
ipp=ipp_p
ipp%pierce(:)%IncVtec=ipp%pierce(:)%IncVtec*Factor
!单位化为米
!在整个ionoGrid数据处理内部，VTEC都化为等效的长度单位
!------------------
! call test
!输出中间结果
IGP%BandNum=11
!总带数
IGP%SysType=2
!1为1期，2为试验，3为全球
IGP%week=week
IGP%sow=sow
! call IGP_delay_give
call IGP_delay_givei
do i=0,10
IGP%GridBand(i)%BandID=i
!带号
do j=1,210
iGP%GridBand(i)%Grid(j)%GridID=j
!格网点编号
end do
end do
!-------
iGP_P=iGP
ipp_p=ipp
!-------
return
end subroutine ionoGrid
subroutine ionogrid(iGP_p,ipp_p,week,sow,Grid_p)
!---------------------------------subroutine comment
! Purpose : 由穿刺点计算格网点延迟及GIVEI
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-12 15:26:21
! 周和秒由外部传入
!
! *** 说明：需要配置文件 gird.dat 不能删除
! 该函数为计算主函数
!-----------------------------------------------------
! Input Parameters :
! ipp-----------------穿刺点结构体
! week----------------格网点周
! sow-----------------格网点秒
! week 和 sow 由外部传入
! Grid_p-----------------格网结构
! Output Parameters :
! IGP-----------------格网点结构体
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use iono_grid
use def_fi
implicit real*8(a-h,o-z)
real*8,parameter:: Factor = 0.1651d0
!基于B1频点 TECU 到 长度单位（米） 转换系数
integer *2 week
integer *4 sow
type (Ipplist) :: iPP_p
type (IGPlist) :: iGP_p
! type (gridpoint) :: Grid_p
! type (bandsgrid) :: grid_p(0:10)
type (bandsgrid) :: grid_p
!type (gridpoint):: grid_song(202,0:10)
! Grid=Grid_p
!该部分由FORTRAN读入，并放在公共内存，这里传回给Grid
!----------------------------------------------------
do i=0,10
do j=1,202
Grid(j,i)%xlong=Grid_p%bandpoint(i)%apoint(j)%xlong
Grid(j,i)%xlat =Grid_p%bandpoint(i)%apoint(j)%xlat
end do
end do
!-----------------------------------------------------
!------------------
ipp=ipp_p
ipp%pierce(:)%IncVtec=ipp%pierce(:)%IncVtec*Factor
!单位化为米
!在整个ionoGrid数据处理内部，VTEC都化为等效的长度单位
!------------------
! call test
!输出中间结果
IGP%BandNum=11
!总带数
IGP%SysType=2
!1为1期，2为试验，3为全球
IGP%week=week
IGP%sow=sow
! call IGP_delay_give
call IGP_delay_givei
do i=0,10
IGP%GridBand(i)%BandID=i
!带号
do j=1,210
iGP%GridBand(i)%Grid(j)%GridID=j
!格网点编号
end do
end do
!-------
iGP_P=iGP
ipp_p=ipp
!-------
return
end subroutine ionoGrid
subroutine IGP_delay_givei
!---------------------------------subroutine comment
! Purpose : 计算格网点的delay及giveI
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-14 02:10:42
!
!-----------------------------------------------------
use def_fi
implicit real*8(a-h,o-z)
integer::ippband(ipp%n_pierce+1,0:10)
call ipp2band(ippband,ipp%n_pierce)
!把穿刺点分到各个带中
!ippband记录各个带中穿刺点序号
call IGP_delay(ippband)
!由穿刺点计算格网点延迟
!call IPP_cVtec(ippband)
call IPP_cVTEC
!由格网点反算穿刺点延迟
call IGP_delay1(ippband)
!第二次计算，已经剔除了穿刺点VTEC与CVTEC误差大于 TOL的数据
call IGP_givei(ippband)
!计算GIVEI值
!do i=1,ipp%n_pierce
!
! write(*,*) 'routine from test1',ipp%n_pierce
! write(*,101) ipp%n_pierce, ipp%pierce(i)%sow,ipp%pierce(i)%phi, ipp%pierce(i)%lam ,&
! ipp%pierce(i)%IncVtec,ipp%pierce(i)%cVtec
!
!end do
!
!101 format(2I8,4F16.6)
end subroutine IGP_delay_givei
subroutine IPP_cVTEC
!---------------------------------subroutine comment
! Purpose : 由格网点反算穿刺点理论值（方法2）
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-03 17:30:20
! 本版本为方法2，基本思路是对每一个穿刺点
! 搜索其周围的格网点，并判断格网点delay是否已经解算出
! 当超过三个格网点解算出后，则解算其VTEC理论值
! 该方法比方法1简便
!-----------------------------------------------------
! Input Parameters :
!
!
! Output Parameters :
!
!
! Subroutines used :
! *. 由经纬度搜索格网点函数
! x1=Dmod(-42d0,5d0) !-2
! x2=MODULO(-42d0,5d0) ! 3
! x3=Dmod(42d0,5d0) !2
! x4=modulo(42d0,5d0) !2
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use def_fi
implicit real*8(a-h,o-z)
integer::ippband(ipp%n_pierce+1,0:10)
ipp%pierce(:)%cVtec = 0d0
!设置所有穿刺点理论值初值为0
do i=1,ipp%n_pierce
xlongP=ipp%pierce(i)%lam
xlatP=ipp%pierce(i)%phi
!正负55度以内跨度5度，高纬度跨度10度
if(xlatP<=55.AND.xlatP>=-55) then
span=5d0
else if (xlatP>55.or.xlatP<-55) then
span=10d0
end if
remLong=MODULO(xlongP,span)
remLat=MODULO(xlatP,span)
!MODULO范例见程序开头注释部分
GridLong1=xlongP-remLong
GridLong2=GridLong1+span
GridLat1=xlatP-remlat
GridLat2=GridLat1+span
!化到正负180度之间
Glong1= anpm (Gridlong1)
Glong2= anpm(Gridlong2)
GLat1= anpm(Gridlat1)
Glat2= anpm(Gridlat2)
! grid_search(iflag,iband,iindex,iband2,iindex2,Glong1,Glat1,ierr)
isum=0
vtectmp=0d0
wtmp=0d0
!下面开始搜索格网点编号
call grid_search(iflag,ib1,ig1,iband2,iindex2,Glong1,Glat1,ierr)
if(iflag/=0) then
if (IGP%GridBand(ib1)%mask(ig1)==1) then
call weight(w,Glong1,Glat1,xlongP,xlatP)
!计算格网点与穿刺点之间的权重
!vtectmp=vtectmp+w*ipp%pierce(itmp2)%IncVtec
vtectmp=vtectmp+w*iGP%GridBand(ib1)%Grid(ig1)%delay
!可以再加模型改正部分
wtmp=wtmp+w
isum=isum+1
end if
end if
call grid_search(iflag,ib2,ig2,iband2,iindex2,Glong1,Glat2,ierr)
if(iflag/=0) then
if (IGP%GridBand(ib2)%mask(ig2)==1) then
call weight(w,Glong1,Glat2,xlongP,xlatP)
!计算格网点与穿刺点之间的权重
!vtectmp=vtectmp+w*ipp%pierce(itmp2)%IncVtec
vtectmp=vtectmp+w*iGP%GridBand(ib2)%Grid(ig2)%delay
!可以再加模型改正部分
wtmp=wtmp+w
isum=isum+1
end if
end if
call grid_search(iflag,ib3,ig3,iband2,iindex2,Glong2,Glat1,ierr)
if(iflag/=0) then
if (IGP%GridBand(ib3)%mask(ig3)==1) then
call weight(w,Glong2,Glat1,xlongP,xlatP)
!计算格网点与穿刺点之间的权重
!vtectmp=vtectmp+w*ipp%pierce(itmp2)%IncVtec
vtectmp=vtectmp+w*iGP%GridBand(ib3)%Grid(ig3)%delay
!可以再加模型改正部分
wtmp=wtmp+w
isum=isum+1
end if
end if
call grid_search(iflag,ib4,ig4,iband2,iindex2,Glong2,Glat2,ierr)
if(iflag/=0) then
if (IGP%GridBand(ib4)%mask(ig4)==1) then
call weight(w,Glong2,Glat2,xlongP,xlatP)
!计算格网点与穿刺点之间的权重
!vtectmp=vtectmp+w*ipp%pierce(itmp2)%IncVtec
vtectmp=vtectmp+w*iGP%GridBand(ib4)%Grid(ig4)%delay
!可以再加模型改正部分
wtmp=wtmp+w
isum=isum+1
end if
end if
if (isum>=3) then
ipp%pierce(i)%cVtec=vtectmp/wtmp
else if (isum<3) then
ipp%pierce(i)%cVtec=0d0
end if
end do
end subroutine IPP_cVTEC
subroutine IGP_givei(ippband)
!---------------------------------subroutine comment
! Purpose : calculate the GIVEI of the Grid point
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-02 09:26:35
! calculate the GIVEI
! IT SHOULD BE CALLED AFTHER THE DELAY HAD BEEN CALCULATED
! V1.1------2014-05-05 16:52:03
! 修复未满足反算条件情况下穿刺点统计的BUG
! 并在统计结果上加上中位数统计信息
! V2.0------2014-05-14 15:35:12
! 统计垂直格网点点离散方法
! Ref : Ionospheric Time Delay Estimation using
! IDW Grid Model for GAGAN
! Niranjan Parasad
!-----------------------------------------------------
! Input Parameters :
! ippband--------------
! pierce point's in each band's index
! Output Parameters :
! UPDATE THE GIVEI OF THE GRID POINT WHICH WAS SAVED IN THE MODULE
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use def_fi
implicit real*8(a-h,o-z)
integer::ippband(ipp%n_pierce+1,0:10)
!---------------
!real*8:: Vec(ipp%n_pierce+2)
!最大可能性为所有穿刺点
!一个格网点附近所有穿刺点垂直序列
!增加时间 2014-05-14 15:38:48 宋叶志
!随着版本 V2.0 新引入变量
!integer*8:: indexVec(ipp%n_pierce+2,3)
!垂直误差指标
!第一个指标为 站号
!第二个指标为 卫星号
!第三个指标为 时间序列次数
!---------------
!real*8,parameter:: Factor = 0.1651d0
!基于B1频点 TECU 到 长度单位（米） 转换系数
real*8,parameter::kappa=3.2905d0
!99.9%置信概率分位数
do i=0,10
IGP%GridBand(i)%Grid(:)%GIVEI=15
!设置所有格网点初值为未监测状态
end do
do iband=0,10
itmp=sum(IGP%GridBand(iband)%mask)
!这个带总共解算出多少格网点
if (itmp<5) cycle
! 如果设置 itmp<3 则退出当次循环 ，解释如下
! !如果一个带总解算出的格网点少于3个
! !则这个带不可能由格网点解算出穿刺点延迟
! !那么就退出本带循环
!
! 这里更进一步的这里设置 itmp<5 退出当次循环
! 因为如果解算出的格网点低于5个，则不可能由
! 格网点算出超过3个穿刺点的理论延迟值
! 而少于3个穿刺点的理论延迟，则无法对一个格网点
! 计算其 GIVE值
! 这样的话，计算这个带中穿刺点理论值已经没有必要
! 故而退出本带内穿刺点理论值的计算
do k=1,210
if (IGP%GridBand(iband)%mask(k)==0) cycle
!MASK记录了是否成功解算出delay
!如果格网点未算出 delay则退出 当次循环
!格网点如果已经算出来的话，则其周围至少三个象限内
!都有穿刺点，所以这里不需要再进行是否满足三个象限的判断
!存储 give有关数值
sumtmp1=0d0
sumtmp2=0d0
sumtmp3=0d0
M=0
!记录计算GIVE时用到穿刺点的个数
call grid_index(iband,k,xlongGrid,xlatGrid,ierr)
sum_err=0d0
!本次循环仅仅用于统计误差均值
do k1=1,ipp%n_pierce
if(ipp%pierce(k1)%cVtec==0) cycle
! !如果穿刺点处反算VTEC为0，则不参加解算 deq 2014-05-05 08:41
xlongP=ipp%pierce(k1)%lam
xlatP=ipp%pierce(k1)%phi
!穿刺点经纬度
!----------------------
call quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
! 判断穿刺点是否落在格网点5度范围内的象限
!及其落在第几个象限
if(iq/=0) then
call weight(wei,xlongGrid,xlatGrid,xlongP,xlatP)
!计算权重
err=ipp%pierce(k1)%cVtec-ipp%pierce(k1)%IncVtec
!计算穿刺点理论值与观测值的误差
! sumerr=sumerr+err
!统计误差和
sumtmp1=sumtmp1+err*wei
sumtmp2=sumtmp2+wei
!统计 GIVE 第一项
sumtmp3=sumtmp3+dabs(err)*wei
! M=M+1
end if
end do
!对所有穿刺点k1循环结束
ave_err=sumtmp1/sumtmp2
!实际为M个点满足统计条件
abs_err=sumtmp3/sumtmp2
!统计置信区间 2014-05-13 14:07:31 syz
!==================================================================
!存储 give有关数值
sumtmp1=0d0
sumtmp2=0d0
do k1=1,ipp%n_pierce
if(ipp%pierce(k1)%cVtec==0) cycle
! !如果穿刺点处反算VTEC为0，则不参加解算 deq 2014-05-05 08:41
xlongP=ipp%pierce(k1)%lam
xlatP=ipp%pierce(k1)%phi
!穿刺点经纬度
!----------------------
call quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
! 判断穿刺点是否落在格网点5度范围内的象限
!及其落在第几个象限
if(iq/=0) then
call weight(wei,xlongGrid,xlatGrid,xlongP,xlatP)
!计算权重
err=ipp%pierce(k1)%cVtec-ipp%pierce(k1)%IncVtec
!计算穿刺点理论值与观测值的误差
sumtmp1 = sumtmp1 + (err - ave_err)*(err - ave_err)
sumtmp2 = sumtmp2 + 1d0;
end if
end do ! k1, 对所有穿刺点k1循环结束
!=============================================================
!GIVE=dabs(ave_err)+kappa*dsqrt(sumtmp1/sumtmp2)
GIVE=abs_err+kappa*dsqrt(sumtmp1/sumtmp2)
!GIVE值
call Give2GiveI(iGiveI,Give)
IGP%GridBand(iband)%Grid(k)%GIVEI=iGiveI
end do !k ,对所有格网点循环结束
end do !iband ,对所有带循环结束
end subroutine IGP_givei
subroutine IGP_delay(ippband)
!subroutine IGP_delay(ipp,iGP)
!---------------------------------subroutine comment
! Purpose :
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-06 19:16:33
! 计算所有格网点延迟
!-----------------------------------------------------
! Input Parameters :
! ipp--------穿刺点结构体数组
!
! Output Parameters :
! igp-------格网点数组，其中包含格网点延迟
!
! ippband----穿刺点分布在不同带中的索引
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use def_fi
implicit real*8(a-h,o-z)
integer::ippband(ipp%n_pierce+1,0:10)
!real*8,parameter:: Factor = 0.1651d0
!基于B1频点 TECU 到 长度单位（米） 转换系数
do iband=0,10 !------带循环
!对所有的带循环
!---------------------------------------------
! if (ippband(1,iband)==0) cycle
! !如果该带中第一个指标为0，表示该带就没有穿刺点
! !则退出本次带的所有格网点循环
!---------------------------------------------
!---------------------------------------------
if (ippband(3,iband)==0) cycle
! !如果该带中第3个指标为0，表示该带即便有1个或2个穿刺点也不足以
! 解算出一个格网点
! 这样避免了一些带只有一两个穿刺点，而造成本带进行循环
!---------------------------------------------
inumber=201
!每个带格网点个数
if(iband==8) inumber=200
if(iband==9) inumber=192
if(iband==10) inumber=192
igp%GridBand(iband)%Mask(:)=0
!一个带的格网点掩码初值设置为0
do k1=1,inumber !-----对每一个格网点
call grid_index(iband,k1,xlonggrid,xlatgrid,ierr)
! gridlong ---格网点经度
! gridlat ----格网点纬度
! 根据带号和格网编号，给出格网点的经纬度
it1=0
it2=0
it3=0
it4=0
! 4个象限的循环指标
do k2=1,ipp%n_pierce+1
!------------对所有穿刺点循环
!------------目的是判断每个格网点四个象限内穿刺点的个数
! 虽然表面上N很大，但是每次循环到每个带中无穿刺点时，自动结束
! 所以循环量较小
itmp=ippband(k2,iband)
!该带中实际的穿刺点指标
!索引这个带(第iband带)中的第k2个穿刺点的序号
!从该带中第一个穿刺点索引，索引到为0表示这个带中已经没有穿刺点
if(itmp==0) exit
!碰到第一个为0时候，表示该带中已经没有穿刺点，退出循环
xlongP=ipp%pierce(itmp)%lam
xlatP=ipp%pierce(itmp)%phi
!穿刺点经纬度
call quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
! 判断穿刺点是否落在格网点5度范围内的象限
!及其落在第几个象限
!同一个象限的穿刺点个数不累计，这里仅仅是为了判断
!是否满足格网点多于3个象限
select case (iq)
case (1)
it1=1
case (2)
it2=1
case (3)
it3=1
case (4)
it4=1
end select
end do !-k2 循环结束，一个格网点对带内所有穿刺点循环结束
ittmp=it1+it2+it3+it4
if (ittmp>=3) then
!大于等于3个象限时解算
vtectmp=0d0 !vtec
wtmp=0d0 !权重
do k3=1,ipp%n_pierce+1
itmp2=ippband(k3,iband)
!该带中实际的穿刺点指标
if(itmp2==0) exit
!碰到第一个为0时候，表示该带中已经没有穿刺点，退出循环
xlongP=ipp%pierce(itmp2)%lam
xlatP=ipp%pierce(itmp2)%phi
!穿刺点经纬度
call quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
! 判断穿刺点是否落在格网点5度范围内的象限
!及其落在第几个象限
if(iq/=0) then
call weight(w,xlongGrid,xlatGrid,xlongP,xlatP)
!计算格网点与穿刺点之间的权重
vtectmp=vtectmp+w*ipp%pierce(itmp2)%IncVtec
!可以再加模型改正部分
wtmp=wtmp+w
end if
!X36B_d%GridBand(iband)%IPG(k1)%GridDelayCrt= vtectmp/wtmp
tmp_delay=vtectmp/wtmp
iGP%Gridband(iband)%Grid(k1)%delay=tmp_delay
!延迟值
!iGP%Gridband(iband)%Mask(k1)=1
!如果可以算出delay值，则掩码为1
end do
iGP%Gridband(iband)%Mask(k1)=1
!如果可以算出delay值，则掩码为1
!改到循环体外 2014-06-12
end if !如果象限大于3个的条件及条件内处理结束
end do ! k1 带内格网点循环结束
end do ! iband 所有带循环结束
end subroutine IGP_delay
subroutine IGP_delay1(ippband)
!---------------------------------subroutine comment
! Purpose : 计算格网点延迟，第二次计算
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-06 19:16:33
! 计算所有格网点延迟
! V2.0-----2014-05-29 15:34:06
! 如果 穿刺点 |VTEC-CVTEC|> TOL 则这一类穿刺点
! 第二次计算时候不参与 格网点计算
!-----------------------------------------------------
! Input Parameters :
! ipp--------穿刺点结构体数组
!
! Output Parameters :
! igp-------格网点数组，其中包含格网点延迟
!
! ippband----穿刺点分布在不同带中的索引
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use def_fi
implicit real*8(a-h,o-z)
integer::ippband(ipp%n_pierce+1,0:10)
real*8,parameter :: TOL =5D0
!用于判断穿刺点理论值与观测值误差容限
!real*8,parameter:: Factor = 0.1651d0
!基于B1频点 TECU 到 长度单位（米） 转换系数
do iband=0,10 !------带循环
!对所有的带循环
!---------------------------------------------
! if (ippband(1,iband)==0) cycle
! !如果该带中第一个指标为0，表示该带就没有穿刺点
! !则退出本次带的所有格网点循环
!---------------------------------------------
!---------------------------------------------
if (ippband(3,iband)==0) cycle
! !如果该带中第3个指标为0，表示该带即便有1个或2个穿刺点也不足以
! 解算出一个格网点
! 这样避免了一些带只有一两个穿刺点，而造成本带进行循环
!---------------------------------------------
inumber=201
!每个带格网点个数
if(iband==8) inumber=200
if(iband==9) inumber=192
if(iband==10) inumber=192
igp%GridBand(iband)%Mask(:)=0
!一个带的格网点掩码初值设置为0
do k1=1,inumber !-----对每一个格网点
call grid_index(iband,k1,xlonggrid,xlatgrid,ierr)
! gridlong ---格网点经度
! gridlat ----格网点纬度
! 根据带号和格网编号，给出格网点的经纬度
it1=0
it2=0
it3=0
it4=0
! 4个象限的循环指标
do k2=1,ipp%n_pierce+1
!------------对所有穿刺点循环
!------------目的是判断每个格网点四个象限内穿刺点的个数
! 虽然表面上N很大，但是每次循环到每个带中无穿刺点时，自动结束
! 所以循环量较小
!======================================
if (dabs(ipp%pierce(k2)%cVtec-ipp%pierce(k2)%IncVtec)> TOL) cycle
!如果穿刺点理论值与测量值误差绝对值大于误差剔除标准
!则剔除这一类穿刺点不参与计算
!2014-05-29 15:06:33 song.yz
!======================================
itmp=ippband(k2,iband)
!该带中实际的穿刺点指标
!索引这个带(第iband带)中的第k2个穿刺点的序号
!从该带中第一个穿刺点索引，索引到为0表示这个带中已经没有穿刺点
if(itmp==0) exit
!碰到第一个为0时候，表示该带中已经没有穿刺点，退出循环
xlongP=ipp%pierce(itmp)%lam
xlatP=ipp%pierce(itmp)%phi
!穿刺点经纬度
call quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
! 判断穿刺点是否落在格网点5度范围内的象限
!及其落在第几个象限
!同一个象限的穿刺点个数不累计，这里仅仅是为了判断
!是否满足格网点多于3个象限
select case (iq)
case (1)
it1=1
case (2)
it2=1
case (3)
it3=1
case (4)
it4=1
end select
end do !-k2 循环结束，一个格网点对带内所有穿刺点循环结束
ittmp=it1+it2+it3+it4
if (ittmp>=3) then
!大于等于3个象限时解算
vtectmp=0d0 !vtec
wtmp=0d0 !权重
do k3=1,ipp%n_pierce+1
itmp2=ippband(k3,iband)
!该带中实际的穿刺点指标
if(itmp2==0) exit
!碰到第一个为0时候，表示该带中已经没有穿刺点，退出循环
xlongP=ipp%pierce(itmp2)%lam
xlatP=ipp%pierce(itmp2)%phi
!穿刺点经纬度
call quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
! 判断穿刺点是否落在格网点5度范围内的象限
!及其落在第几个象限
if(iq/=0) then
call weight(w,xlongGrid,xlatGrid,xlongP,xlatP)
!计算格网点与穿刺点之间的权重
vtectmp=vtectmp+w*ipp%pierce(itmp2)%IncVtec
!可以再加模型改正部分
wtmp=wtmp+w
end if
!X36B_d%GridBand(iband)%IPG(k1)%GridDelayCrt= vtectmp/wtmp
tmp_delay=vtectmp/wtmp
iGP%Gridband(iband)%Grid(k1)%delay=tmp_delay
!延迟值
!iGP%Gridband(iband)%Mask(k1)=1
!如果可以算出delay值，则掩码为1
end do
iGP%Gridband(iband)%Mask(k1)=1
!如果可以算出delay值，则掩码为1
!改到循环体外 ，2014-06-12
end if !如果象限大于3个的条件及条件内处理结束
end do ! k1 带内格网点循环结束
end do ! iband 所有带循环结束
end subroutine IGP_delay1
subroutine quadrant(iq,xlongGrid,xlatgrid,xlongP,xlatP)
!---------------------------------subroutine comment
! Purpose : 根据格网点经纬度和穿刺点经纬度，判断穿刺点位于格网点的
! 的哪各象限
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-06 23:14:58
! 根据格网点经纬度和穿刺点经纬度
! 判断穿刺点位于格网点的哪个象限
! 特别需要注意的是在格网点经度是+/180时的特殊边界情况
!
! 一种看似可行的方法，其实不可行，即对格网点和穿刺点的经度
! 同时加 360的倍数 如720
! 这样同样要面临，起始点的问题
! 应该按照本程序中的方法进行处理
! V2.0 ----2014-06-12
! 对于落在象限的经纬度范围由原先的5度改为可设置
! 这里设置为2.5度
!-----------------------------------------------------
! Input Parameters :
! xlonggrid-----格网点经度
! xlatgrid ------格网点纬度
! xlongP-------穿刺点经度
! xlongP-------穿刺点纬度
! Output Parameters :
! iq -- 0----穿刺点不在格网点5度附近
! 1--- 穿刺点位于格网点第一象限 [0，90]
! 2----穿刺点位于格网点第二象限 [90，180]
! 3----穿刺点位于格网点第三象限 [180，270]
! 4----穿刺点位于格网点第四象限 [270，360]
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8,parameter:: xp = 5d0 !正数
real*8,parameter:: xm = -5d0 !负数
!-----------非常重要
if (xlongGrid==-180.and.(xlongP>175.and.xlongP<180)) then
xlongP=xlongp-360
end if
!处理经度在 -180度边界线上的交叉问题
if (xlongGrid==180.and.(xlongP>-180.and.xlongP<-175)) then
xlongP=xlongP+360
end if
!处理经度在 180度边界线上的交叉问题
tmp1=xlongP-xlongGrid
tmp2=xlatP-xlatGrid
iq=0
!设置初状态
if ((tmp10d0).and.(tmp20d0)) then
iq=1
!第一象限
else if ((tmp10d0).and.(tmp2<0d0.and.tmp2>xm)) then
iq=4
else if ((tmp1>=xm.and.tmp1<=0d0) .and. (tmp2<=xp.and.tmp2>=0d0)) then
iq=2
else if ((tmp1>=xm.and.tmp1<=0d0) .and. (tmp2<=0d0.and.tmp2>=xm)) then
iq=3
end if
end subroutine quadrant
subroutine quadrant1(iq,xlongGrid,xlatgrid,xlongP,xlatP)
!---------------------------------subroutine comment
! Purpose : 根据格网点经纬度和穿刺点经纬度，判断穿刺点位于格网点的
! 的哪各象限
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-06 23:14:58
! 根据格网点经纬度和穿刺点经纬度
! 判断穿刺点位于格网点的哪个象限
! 特别需要注意的是在格网点经度是+/180时的特殊边界情况
!
! 一种看似可行的方法，其实不可行，即对格网点和穿刺点的经度
! 同时加 360的倍数 如720
! 这样同样要面临，起始点的问题
! 应该按照本程序中的方法进行处理
!-----------------------------------------------------
! Input Parameters :
! xlonggrid-----格网点经度
! xlatgrid ------格网点纬度
! xlongP-------穿刺点经度
! xlongP-------穿刺点纬度
! Output Parameters :
! iq -- 0----穿刺点不在格网点5度附近
! 1--- 穿刺点位于格网点第一象限 [0，90]
! 2----穿刺点位于格网点第二象限 [90，180]
! 3----穿刺点位于格网点第三象限 [180，270]
! 4----穿刺点位于格网点第四象限 [270，360]
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
implicit real*8(a-h,o-z)
!-----------非常重要
if (xlongGrid==-180.and.(xlongP>175.and.xlongP<180)) then
xlongP=xlongp-360
end if
!处理经度在 -180度边界线上的交叉问题
if (xlongGrid==180.and.(xlongP>-180.and.xlongP<-175)) then
xlongP=xlongP+360
end if
!处理经度在 180度边界线上的交叉问题
tmp1=xlongP-xlongGrid
tmp2=xlatP-xlatGrid
iq=0
!设置初状态
if ((tmp1<5d0.and.tmp1>0d0).and.(tmp2<5d0.and.tmp2>0d0)) then
iq=1
!第一象限
else if ((tmp1<5d0.and.tmp1>0d0).and.(tmp2<0d0.and.tmp2>-5d0)) then
iq=4
else if ((tmp1>=-5d0.and.tmp1<=0d0) .and. (tmp2<=5d0.and.tmp2>=0d0)) then
iq=2
else if ((tmp1>=-5d0.and.tmp1<=0d0) .and. (tmp2<=0d0.and.tmp2>=-5d0)) then
iq=3
end if
end subroutine quadrant1
subroutine weight(w,xlongGrid,xlatGrid,xlongP,xlatP)
!---------------------------------subroutine comment
! Purpose : 计算格网点与穿刺点的权重
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-07 00:00:08
! 根据球面距离计算权重
! 注意输入参数均为角度（算权重时要化为弧度）
! V2.0----2014-06-12 18:49:12
! 去掉球半径 6378.1改用单位1， 半径作为常数，用在不同的权上，最终消去了。
!-----------------------------------------------------
! Input Parameters :
! xlongGrid------格网点经度（角度）
! xlatGrid-------格网点纬度（角度）
! xlongP--------穿刺点经度（角度）
! xlatP---------穿刺点纬度（角度）
! Output Parameters :
! w------------权重
!
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
implicit real*8(a-h,o-z)
!tmp = dsin(xlatGrid)*dsin(xlatP)+dcos(xlatGrid)*dcos(xlatP)*dcos(xlongGrid-xlongP)
arc=0.01745329251994d0
!degree two arc ,arc = pi/180
tmp = dsin(xlatGrid*arc)*dsin(xlatP*arc)+dcos(xlatGrid*arc)*dcos(xlatP*arc)*dcos(xlongGrid*arc-xlongP*arc)
!d = 6378.1*dacos(tmp)
!w = 1d0/d
eps=3d-4
!加一小常数，防止奇异
d= dacos(tmp)+eps
w=1d0/d
end subroutine weight
subroutine ipp2band(ippband,N)
!---------------------------------subroutine comment
! Purpose : set the piece points to different bands
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-31 16:20:04 （had been tested）
! set the piece points to different bands
! 用ippband记录各个带中穿刺点的序列
! 如第0个带穿刺点序列为 7，15，24，...
! ...
! 第10带穿刺点序列为 8，12，36 ...
!
! 经过如此存储后，允许同一个穿刺点在不同的带中被记录
! 如被经度交叉的带中同时记录
!-----------------------------------------------------
! Input Parameters :
! ipp --------piece points
! N-----------the number of iono piece points
! Output Parameters :
! ippband-------piece points in different bands
!
! 0 1 2 3 4 .... 10
! . . . . . .
! . . . . . .
!
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use def_fi
implicit real*8(a-h,o-z)
!include 'def.fi'
!!引入数据结构
!type (ipplist) :: ipp
integer::ippband(N+1,0:10)
integer::k(0:10)
! the index of every band
! do i=1,5
!
! write(*,*) 'routine from ipp2band',ipp%n_pierce
! write(*,101) ipp%n_pierce, ipp%pierce(i)%sow,ipp%pierce(i)%phi, ipp%pierce(i)%lam ,&
! ipp%pierce(i)%IncVtec,ipp%pierce(i)%cVtec
!
!end do
!
!101 format(2I8,4F16.6)
ippband=0
k(0:10)=1
do i=1,ipp%n_pierce
xlong=ipp%pierce(i)%lam
xlat=ipp%pierce(i)%phi
!longitude and latitude of pierce point
xlong=anpm(xlong)
!transfer the longitude to 0 - 180(degree)
!-----------------------------------------
do iband=0,8
! !cross 5 degree in both plus and minus direction **** important
tmp1=-180+40*iband-5
tmp2=tmp1+40+10
if (xlong>=tmp1.and.xlong175) then
ippband(k(0),0)=i
k(0)=k(0)+1
end if
! 175-----180 也化到第0带
if (xlong>-180.and.xlong<-175) then
ippband(k,8)=i
k(8)=k(8)+1
end if
! -180 ---- -175 也化到第8带
!-----------------------------------------
if(xlat>55) then
!向下交叉5度
ippband(k,9)=i
k(9)=k(9)+1
end if
if(xlat<-55) then
!向上交叉5度
ippband(k,10)=i
k(10)=k(10)+1
end if
end do
end subroutine ipp2band
subroutine grid_index(iband,iindex,xlong,xlat,ierr)
!---------------------------------subroutine comment
! Purpose : Output the longitude and latitude by the band ID and the point ID
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-01 16:51:06
! Output the longitude and latitude by the band ID and the point ID
! V2.0-----2014-03-04 18:01:22
! get the longitude and latitude directly by new data structure
!-----------------------------------------------------
! Input Parameters :
! iband-----------band ID
! iindex----------point ID (in the band)
! Output Parameters :
! xlong-----------longitude
! xlat------------latitude
! ierr------------flag for record whether it is wrong
! it means the wrong ID if ierr equals minus one
! Subroutines used :
! *.
! *.
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use iono_grid
implicit real*8(a-h,o-z)
if(iindex<1.OR.iindex>201) ierr=-1
if(iband>10.OR.iband<0) ierr=-1
if(iband==8.AND.iindex>200) ierr=-1
!there are 200 Grid points in band 8
xlong=grid(iindex,iband)%xlong
xlat=grid(iindex,iband)%xlat
!output the information of the Grid
! do i=0,2
! do j=1,202
! write(*,*) Grid(j,i)%xlong, Grid(j,i)%xlat,'Grid point'
! end do
! end do
!-----------------------------------------------------
end subroutine grid_index
subroutine grid_search(iflag,iband,iindex,iband2,iindex2,xlong0,xlat0,ierr)
!---------------------------------subroutine comment
! Purpose : seach the Grid ID by longitude and latitude
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-03-01 10:46:21
! 根据经纬度搜索格网点
! 是否搜索到以 is,is2作为判断标准
! V2.0 ---合并 is与 is2 ，以 iflag作为输出指标
! iflag 意义见 output parameters
! V3.0 ---2014-03-04 17:23:20
! 采用格网点结构体二维矩阵索引，这样就不需要计算每个格网点在所有矩阵中的绝对指标
! 相应的模块也作升级为 iono_grid
! V3.1 ---2014-03-05 17:19:59
! 判断两个角度是否相等，不以 绝对相等为依据
! 而是当两者之差绝对值小于 一个很小的数 eps(设置为1D-15)时，认为两者相等
! 这样避免了，角度和弧度相互转化之后，带来的精度影响最终的匹配
! V3.2 ---2014-03-06 11:29:32
! Normalize angle into the range -180 <= A < +180
! 不再度输入角度检查是否在 +/-180之间
! 而是把输入的经纬度化到 +/-180
! 如果纬度 在 +/- 90之外，则 令 ierr为 -1
!-----------------------------------------------------
! Input Parameters :
! xlong---------经度（单位：角度） 范围 -180---180
! xlat----------纬度（单位：角度） 范围 -90----90
!
! Output Parameters :
! iflag 0 未匹配到格网点
! 1 匹配到 1-8带中某格网点，但未匹配到9或10带
! 带号及带内号见 iband, iindex
! 2 匹配到 9 或 10带，但未匹配到1-8带
! 带号及带内号见 iband2,iindex2
! 3 同时匹配到1-8和9（或10）
! 带号和带内号分别见
! iband,iindex,iband2,iindex2
! ==================================================
! iband--------如果落在1-8某带中，则输出带编号
! 无论是否匹配到穿刺点，iband一定会输出
! 而是否匹配成功，则以 is 为依据
! 因为角度总在-180--180之间，所以一定iband一定有值
! iindex-------输出该iband带内号
! 如果匹配到穿刺点，则输出该编号
! 如果未匹配到穿刺点，则输出该带的最后一个编号
! 未匹配到，则输出-1
! iband2-------如果既落在1-8某带中，也落在9-10某带中，则输出该带号
! (9或10)
! 如果未落在 9 或者 10带中
! iband2=-1
! iindex2------该带 iband2（9或10） 的带内编号
! 未匹配到输出-1
! ierr---------ierr为-1，则输入数据范围不在区间内，不予计算
! Subroutines used :
! *. band_interval(iband,ibegin,iend,ierr) 输出各带的始末指标
! *.
!-----------------------------------------------------
! Notes :
! 1.无论是否匹配到穿刺点
! 任何一个点的经纬度范围必然落在1-8带中，所以1-8带 iband一定会输出
! 如果一个点经纬度范围也落在9或10带中，则iband2会输出9或10
!
! 2. 如果为落到9或10区域，则 iband2=-1, iindex2=202
! 这里 iband2默认不取0，是因为有第0带，防止引起歧义
!
! 3. 对于1-8区域，如果未匹配上，则 index=-1
!
! 4. 对于输出结果首先要看 iflag 来判断匹配结果，然后再看其他量
!
! Ref. :
! MOPS for GPS WAAS DO-229D 2006 P252
!
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
use iono_grid !gridMAT
implicit real*8(a-h,o-z)
real*8,parameter::eps=1d-15
!--------------------------
iband=-1
iindex=-1
iband2=-1
iindex2=-1
!--------------------------
!
iflag=0
is=0
is2=0
!-----------------------
! 先把经纬度化为 +/- 180之间
! 如果纬度 在 +/- 90之外 则不计算，并返回 ierr=-1 作为警告信息
xlong=anpm(xlong0)
xlat=anpm(xlat0)
if (xlat>90.or.xlat<-90) then
ierr=-1
return
end if
!=======================================================================
!=======================================================================
! check the point in which band (1-8) by longitude
do i=0,8
tmp1=-180+40*i
tmp2=tmp1+40
if (xlong>=tmp1.and.xlong=60) then
iband2=9
!第二个带指标为9
! iindex2=0
do i=1,192
! if (xlong==grid(i,9)%xlong.and.xlat==grid(i,9)%xlat) then
! V 3.0
if (dabs(xlong-grid(i,9)%xlong) < eps .and. dabs(xlat-grid(i,9)%xlat) < eps) then
is2=1 !如果和9带中某穿刺点匹配上 is2为1
iindex2=i !
exit !1个带中没有重复，如果搜到格网点则退出该带搜索
end if
end do
end if
if (xlat<=-60) then
iband2=10
!第二个带指标为10
! iindex2=0
do i=1,192
if (dabs( xlong-grid(i,10)%xlong ) < eps .and. dabs(xlat-grid(i,10)%xlat) < eps) then
is2=1 !如果和9带中某穿刺点匹配上 is2为1
iindex2=i !
exit !1个带中没有重复，如果搜到格网点则退出该带搜索
end if
end do
end if ! seach in band 10
! iflag 0 have not get the grid point
! 1 get the grid point in band 1 to 8, but not in
! band 9 or 10. the band ID and the point id
! were recorded by iband and iindex
! 2 the point was seached in band 9 or band 10,but
! not in band 1 to 8.
! the band ID and the point ID were recorded by
! iband2 and iindex2.
! 3 the point was in band 1 to 8 and 9(or 10)
! the band ID and the point ID were recorded by
! iband,iindex,iband2,iindex2
if (is==0.and.is2==0) then
iflag=0
return
else if (is==1.and.is2==0) then
iflag=1
return
else if (is==0.and.is2==1) then
iflag=2
return
else if (is==1.and.is2==1) then
iflag=3
end if
end subroutine grid_search
double precision function anpm ( a )
!---------------------------------subroutine comment
! purpose : Normalize angle into the range -180 <= A < +180
! author : song yezhi
! versions and changes :
! v1.0 ----2014-03-04 18:17:59
! if need to normalize the angle to -pi ---pi
! could set
! dpi = 3.1415926535897932384626d0
! d2pi = 6.283185307179586476925287d0
! Ref . sofa iau_anpm
! v2.0 ---
! if the angle equals 180 after it
! had been normalized to +/-180, set it to -180
!-----------------------------------------------------
! input parameters :
! a-----input angle (unit: degree)
! output parameters :
! anpm-----
! angle in range +/-180
! subroutines used :
! *.
! *.
!-----------------------------------------------------
! copyrigt (c) : (all rights reserved.)
! center for astro-geodynamics
! shanghai astronomical observatory
! chinese academy of sciences , 2014
!----------------------------------------------------
implicit none
double precision a
double precision dpi
parameter ( dpi = 180 )
double precision d2pi
parameter ( d2pi = 360 )
double precision w
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
w = mod(a,d2pi)
if ( abs(w) .ge. dpi ) w = w - sign(d2pi,a)
anpm = w
!
if (anpm==180) anpm=-180
!
end function anpm
subroutine Give2GiveI(iGiveI,Give)
!---------------------------------subroutine comment
! Purpose : Give to GiveI EVALUATION OF GIVEI
! Author : Song Yezhi
! Versions and Changes :
! V1.0 ----2014-04-02 13:41:59
! GIVE to GiveI
! Ref : MOPS for GPS Wass DO-229D 2006
! Appendix A Page A-33
!-----------------------------------------------------
! Input Parameters :
! GIVE : Meters
!
! Output Parameters :
! GIVEI : 0--15 ranks
!
! sigma^2
! GIVEI GIVE Meters GIVE Meters^2
! 0 0.3 0.0084
! 1 0.6 0.0333
! 2 0.9 0.0749
! 3 1.20 0.1331
! 4 1.5 0.2079
! 5 1.8 0.2994
! 6 2.1 0.4075
! 7 2.4 0.5322
! 8 2.7 0.6735
! 9 3.0 0.8315
! 10 3.6 1.1974
! 11 4.5 1.8709
! 12 6.0 3.3260
! 13 15.0 20.7870
! 14 45.0 187.0826
! 15 Not Monitored Not Monitored
!-----------------------------------------------------
! Copyrigt (C) : (All rights reserved.)
! Center for Astro-geodynamics
! Shanghai Astronomical Observatory
! Chinese Academy of Sciences , 2014
!----------------------------------------------------
implicit real*8(a-h,o-z)
real*8::rank(0:14)
rank=(/0.3d0,0.6d0,0.9d0,1.2d0,1.5d0,1.8d0, &
2.1d0,2.4d0,2.7d0,3.0d0,3.6d0,4.5d0, &
6.0d0,15.d0,45.d0 /)
if (Give<=rank(0)) iGiveI=0
do i=1,14
if( give<=rank(i).and.give>rank(i-1)) iGiveI=i
end do
if (Give>=45) IGiveI=15
end subroutine Give2GiveI