请帮忙查看程序，如何用GRads画规则圆和不规则圆

雾由心生

我发现程序最关键的问题是在第二个大循环的内循环中，在确定了经纬度时间层次此四维后不能有效提取到实时的att用于下面紧接着的两个公式的计算，哪位大写能帮忙解决，非常感谢

'reinit'

'open f:/fog/model.ctl'

*  PI_180: PI/180

*  olon: original longitude

*  olat: original latitude

*  bear: bearing

*  dist: distant

*  deltd: delta distance,step length

*  tdis: target distance

*  tlon: target longitude

*  tlat: target latitude

*  att: attenuation

'define PI180=0.017453292519943'

'define olon=122.5'

'define olat=32'

'define tlat=olat'

'define tlon=olon'

'define att=25*qcwsfc'

'define angle=0'

'define fu=0'

'define integ=0'

'define fid=0'

'set font 1'

'set mpdset hires'

'set mproj scaled'

'set gxout shaded'

'set grads off'

'set grid off'

'set t 13'

'set lat 31.7 32.3'

'set lon 122 123'

'd qcwsfc'

'cbarn'

'q w2xy 122.5 32'

x =subwrd (result,3)

y=subwrd (result,6)

'draw mark 3 'x' 'y' 0.1'

*Draw the 25km range circle

bear=0

tdis=25

while (bear<360)

  'angle='bear'*'PI180''

  'fu='tdis'*cos(angle)/60.0/1.86'

  'tlat='olat'+fu'

  'tlon='olon'+'tdis'*sin(angle)/60.0/cos((tlat+fu/2)*PI180)'

  'q define '

  tlat=sublin(result,9)

  tlat=subwrd (tlat,2)

  tlon=sublin(result,10)

  tlon=subwrd (tlon,2)

  'q w2xy 'tlon' 'tlat''

  x=subwrd (result,3)

  y=subwrd (result,6)

  'draw mark 3 'x' 'y' 0.05'

  bear=bear+1

endwhile

*Draw the valid detective distance circle

bear=0

while ( bear<360)

  ' integ=0'

  ' fid=0'

  'dist=25'

  deltd=0.1

  dist=20

  tdis=deltd

*This internal circulation is used to get the tlat and tlon when tdis equals dist

*Get the longitude and latitude of the valid detective distance point

  while(tdis<=dist)

    'angle='bear'*'PI180''

    'fu='tdis'*cos(angle)/60/1.86'

    'tlat='olat'+fu'

    'tlon='olon'+'tdis'*sin(angle)/60/cos((tlat+fu/2)*'PI180')'

*Get the attenuation at tlat and tlon,then compute the integral and fi(d)

    'q define '

    tlat=sublin(result,8)

    tlat=subwrd (tlat,2)

    tlon=sublin(result,9)

    tlon=subwrd (tlon,2)

    'set lon ' tlon

    'set lat ' tlat

    'integ=integ+att*'deltd''

    'fid='tdis'*pow(2.718281828,0.115*integ)'

   

    tdis=tdis+'deltd'

  endwhile

*Draw the Dmax point at the tlat and tlon when tdis equals dist through different att

  'q w2xy 'tlon' 'tlat''

  x=subwrd (result,3)

  y=subwrd (result,6)

  'draw mark 3 'x' 'y' 0.15'

  bear=bear+1

endwhile

;

river

没太看懂你的问题是什么意思，但是我感觉总的来说就是循环把经纬度转换成虚页坐标，然后利用draw mark命令画标记。思路不错，我觉得有问题的是赋值语句的单引号运用不对，比如
'angle='bear'*'PI180''
'fu='tdis'*cos(angle)/60.0/1.86'
'tlat='olat'+fu'
  'tlon='olon'+'tdis'*sin(angle)/60.0/cos((tlat+fu/2)*PI180)'
应该是
angle=bear*PI180
fu=tdis*cos('angle')/60.0/1.86
tlat=olat+'fu'
  tlon=olon+tdis*sin('angle')/60.0/cos(('tlat+fu/2')*PI180)
剩下的以此改一下

mofangbao

虽然代码贴的有点长，不过问题挺有意思

雾由心生

river 发表于 2013-5-9 22:03

                               
登录/注册后可看大图

没太看懂你的问题是什么意思，但是我感觉总的来说就是循环把经纬度转换成虚页坐标，然后利用draw mark命令画 ...

谢谢，你说的方法我试过了，不行的啊，应该不是这个问题，如果程序中去掉第二个while大循环其他不变再运行，则可以画出那个25km的圆，如果按你所说去掉循环中的单引号的话会出更多的错，循环会找不到bear、tdis，相反循环中不对fu、angle等用单引号。这个程序随便找一个可以画set gxout shaded的数据都可以运行，只要对里面的经纬度、时间进行设置，然后人依照各变量取代qcwsfc就可以了，感兴趣的话可以试试

雾由心生

　　　我发现最关键的问题是在第二个大循环的内循环中，在确定了经纬度时间层次此四维后不能有效提取到实时的att用于下面紧接着的两个公式的计算，这两个公式应该不是问题所在，难道是因为底图变量是qcwsfc所以提取不出变量att?还有是不是每次循环中q define所列的定义变量顺序一直在变？前后两个大循环中应用sublin函数提取的都是同两个变量，但sublin(result,num)中num可能就不一样，从整个程序运行来看q define所列的定义变量tlat tlon顺序为8,、9，如果整个程序没有后一个大循环其他一切不变，其顺序则为9、10，然而完整程序运行时第一个循环内sublin(result,num)中num必须设置为9、10才能将外圆画出，设置为8、9则不行，不知为什么

雾由心生

不好意思第一次发帖，是有点长了，下次一定注意。
这个程序的关键不是用什么数据，所以如果各位大侠如果感兴趣，自己任意找个数据将变量qcwsfc换掉再将程序里的基本设置更改一下就可以运行了。

mofangbao

不规则圆是什么样子的,是椭圆吗你的问题描述的不太清楚，这个程序正确的话应该是什么样子的 虽然我运行了一下，提示有错误，但是也能出来一个椭圆，所以想知道正确的应该是什么样子，这样才能有个调试的目标。

雾由心生

mofangbao 发表于 2013-5-10 12:02

                               
登录/注册后可看大图

不规则圆是什么样子的,是椭圆吗你的问题描述的不太清楚，这个程序正确的话应该是什么样子的 虽然我运行了一 ...

我对程序进行仔细的推敲，有对第二个while循环的内循环添加了控制语句if  endif,其他地方也作了一些修改，修改后的程序如下：
'reinit'
'clear '
'open f:/fog/model/model.ctl'

'define PI180=0.017453292519943'
'define olon=122.5'
'define olat=32'
'define tlat=olat'
'define tlon=olon'
'define att=25*qcwsfc'
'define angle=0'
'define fu=0'
'define integ=0'
'define fid=0'
'set mpdset hires'
'set mproj scaled'
'set gxout shaded'
'set t 13'
'set lat 31.7 32.3'
'set lon 122.2 122.8'
'd qcwsfc'
'cbarn'
'q w2xy 122.5 32'
x =subwrd (result,3)
y=subwrd (result,6)
'draw mark 3 'x' 'y' 0.1'
*Draw the outside circle
bear=0
tdis=16
while (bear<360)
  'angle='bear'*'PI180''
  'fu='tdis'*cos(angle)/60.0/1.86'
  'tlat=olat+fu'
  'tlon=olon+'tdis'*sin(angle)/60.0/cos((tlat+fu/2)*PI180))'
  'q define '
  tlat=sublin(result,9)
  tlat=subwrd (tlat,2)
  tlon=sublin(result,10)
  tlon=subwrd (tlon,2)
  'q w2xy 'tlon' 'tlat
  x=subwrd (result,3)
  y=subwrd (result,6)
  'draw mark 3 'x' 'y' 0.05'
  bear=bear+1
endwhile
*Draw the inside circle
bear=0
while(bear<=360)
integ=0
fid=0
dist=16
deltd=0.1
tdis=deltd
while(tdis<=16)
    'angle='bear'*PI180'
    'fu='tdis'*cos(angle)/60/1.86'
    'tlat=olat+fu'
    'tlon=olon+'tdis'*sin(angle)/60/cos((tlat+fu/2)*PI180)'
    'q define '
    tlat=sublin(result,9)
    tlat=subwrd (tlat,2)
    tlon=sublin(result,10)
    tlon=subwrd (tlon,2)
    'set lon ' tlon
    'set lat ' tlat
    'integ=integ+(att-10*log10(0.92)/16)*0.1'
    'fid='tdis'*pow(10,0.05*integ)'
    'q define'
    fid=sublin(result,10)
    fid=subwrd (fid,2)
    if(fid>16.6812)
      break
    endif
    tdis=tdis+0.1
  endwhile
  'q w2xy 'tlon' 'tlat''
  x=subwrd (result,3)
  y=subwrd (result,6)
  'draw mark 3 'x' 'y' 0.05'
bear=bear+1
endwhile
;
程序修改后刻画出的图形如下：

f:/fog/model/a.jpg

雾由心生

f:/fog/model/a.jpg

雾由心生

图片在附件中，修改过的程序画出的是图片a，但是应该画出类似于图片b