Embedded Files
下図のように斜めの領域も指定可能。領域平均値を求めたりする際に使用する。
プログラム一式の圧縮ファイル
コンパイル
$ cd src; make; cd ..
設定ファイルの編集
下記の2つのファイルを適宜編集
para.txt
area_test.txt
実行
$ set_domain
Fortranプログラム
set_domain.f90
program set_domain
! Description:
!
! Author: am
!
! Host: aofd30
! Directory: /work3/am/ASCAT/area_ave/set_domain/src
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 09:53 on 10-08-2012.
!
! Reference
! 点が長方形の内側にあるかチェックする
! http://endoyuta.blogspot.jp/2011/09/blog-post_06.html (See below)
! use
! implicit none
character(len=100)::infle_area*100
real xa(4),ya(4)
real,allocatable::xd(:),yd(:)
integer,allocatable::iar(:),jar(:)
real xp,yp
namelist/domain/xw,xe,ys,yn,dx,dy
namelist/area/infle_area
write(*,'(a)')'Program set_domain starts.'
write(*,*)''
open(1,file="para.txt",action="read")
read(1,domain)
read(1,area)
close(1)
nx=int((xe-xw)/dx)+1
ny=int((xe-xw)/dx)+1
allocate(xd(nx),yd(ny))
open(10,file=infle_area,action="read")
do i=1,4
read(10,*)xa(i),ya(i)
enddo
close(10)
print *,"Area:"
print '(A)',infle_area(1:lnblnk(infle_area))
print *,xa(1),xa(2)
print *,ya(1),ya(3)
print *
call count_area(xa,ya,xw,ys,nx,ny,dx,dy,nar)
print '(A,i5)','nar= ',nar
call set_coord(xw,ys,nx,ny,dx,dy,xd,yd)
allocate(iar(nar),jar(nar))
call set_grid_point(xa,ya,xw,ys,nx,ny,dx,dy,nar,iar,jar)
call test_output_grid_point(nx,ny,xd,yd,nar,iar,jar)
write(*,'(a)')'Done program set_domain.'
write(*,*)
end program set_domain
subroutine test_output_grid_point(nx,ny,xd,yd,nar,iar,jar)
integer,intent(in)::nx,ny
real,intent(in)::xd(nx),yd(ny)
integer,intent(in)::nar,iar(nar),jar(nar)
character(len=500)::out_file
out_file="output/test_output_grid.asc"
print '(A,A)',"Output: ",out_file(1:lnblnk(out_file))
open(20,file=out_file)
print *
do n=1,nar
! print *,iar(n),jar(n),xd(iar(n)),yd(jar(n))
write(20,'(2f15.7)')xd(iar(n)),yd(jar(n))
enddo !n
end subroutine test_output_grid_point
subroutine set_coord(xw,ys,nx,ny,dx,dy,xd,yd)
real,intent(in)::xw,ys,dx,dy
integer,intent(in)::nx,ny
real,intent(out)::xd(nx),yd(ny)
do i=1,nx
xd(i)=xw+float(i-1)*dx
enddo !i
do j=1,ny
yd(j)=ys+float(j-1)*dy
enddo !j
end subroutine set_coord
subroutine count_area(xa,ya,xw,ys,nx,ny,dx,dy,nar)
real,intent(in)::xa(4),ya(4),xw,ys
integer,intent(in)::nx,ny
real,intent(in)::dx,dy
integer,intent(out)::nar
k=0
do j=1,ny
do i=1,nx
xp=xw+float(i-1)*dx
yp=ys+float(j-1)*dy
! print "(2f5.1)", xp, yp
op=0.
flagop=0.
do n=1,3
ax=xa(n+1)-xa(n)
ay=ya(n+1)-ya(n)
bx=xp-xa(n)
by=yp-ya(n)
op=outer_product(ax,ay,bx,by)
if(op < 0.0)flagop=-1.
! print "(4f6.2,f10.1)",ax,ay,bx,by,op
enddo
ax=xa(1)-xa(4)
ay=ya(1)-ya(4)
bx=xp-xa(4)
by=yp-ya(4)
op=outer_product(ax,ay,bx,by)
! print "(4f6.2,f10.1)",ax,ay,bx,by,op
if(op < 0.0)flagop=-1.
if(flagop /= -1.)then
k=k+1
endif
enddo !i
enddo !j
nar=k
return
end
subroutine set_grid_point(xa,ya,xw,ys,nx,ny,dx,dy,nar,iar,jar)
real,intent(in)::xa(4),ya(4),xw,ys
integer,intent(in)::nx,ny
real,intent(in)::dx,dy
integer,intent(in)::nar
integer,intent(inout)::iar(nar),jar(nar)
k=0
do j=1,ny
do i=1,nx
xp=xw+float(i-1)*dx
yp=ys+float(j-1)*dy
! print "(2f5.1)", xp, yp
op=0.
flagop=0.
do n=1,3
ax=xa(n+1)-xa(n)
ay=ya(n+1)-ya(n)
bx=xp-xa(n)
by=yp-ya(n)
op=outer_product(ax,ay,bx,by)
if(op < 0.0)flagop=-1.
! print "(4f6.2,f10.1)",ax,ay,bx,by,op
enddo
ax=xa(1)-xa(4)
ay=ya(1)-ya(4)
bx=xp-xa(4)
by=yp-ya(4)
op=outer_product(ax,ay,bx,by)
! print "(4f6.2,f10.1)",ax,ay,bx,by,op
if(op < 0.0)flagop=-1.
if(flagop /= -1.)then
k=k+1
iar(k)=i
jar(k)=j
! xp_test=xw+float(iar(k)-1)*dx
! yp_test=ys+float(jar(k)-1)*dy
! print *,xp_test,yp_test
endif
enddo !i
enddo !j
return
end
!
! Description:
! z-comp of an outer product of vectors, A and B
!
function outer_product(ax,ay,bx,by)
outer_product=ax*by-bx*ay;
end function outer_product
!
! 点が長方形の内側にあるかチェックする
! http://endoyuta.blogspot.jp/2011/09/blog-post_06.html
!
!点が長方形の内側にあるか否かをチェックするには、外積を使えばよい。
!
!ある長方形(10,10),(30,10),(30,30),(10,30)と、ある点(20,20)があるとしよう。このある点がある長方形の内側にあるのか外側にあるのかを判定したい。このような長方形であれば、外積を使うまでもなくチェックが可能だが、長方形が回転している場合や、長方形以外の多角形の場合などになると、簡単ではなくなる。しかし、外積を使えば一発で分かる。外積はこのようにベクトルAとベクトルBの外積を求めた場合、その値がマイナスであればベクトルBはベクトルAの右側にあることになるので、長方形の各辺とある点が全て右側にあることを外積によって確かめればよい。(HTML5のcanvasの場合、座標は左上が原点となるので注意が必要じゃ。一般的な一番左下を原点とする形態であれば、長方形の各辺を時計回りにベクトルにしていけば、点がベクトルの右側にある場合、外積の数値はマイナスになるが、canvasの場合は反時計回りにベクトルにしていくとよい。マイナスになればベクトルの左側にあることにあり、（全ての辺に対して左側にあるのであれば）それはすなわち内側にあるということだ。)
!
!まず、長方形の(30,10)-(10,10)の辺と点の位置関係を確認しよう。長方形の辺をベクトルにすると(-20,0)になる。長方形の辺の始点と点をベクトルにすると、(-10,10)となる。この２つのベクトルの外積を求めると、下記になる。
!-20*10-10*0 = -200
!
!マイナスなので点は辺の左側にあるということだ。これを全ての辺でチェックして、全ての結果がマイナスであれば内側にあるということになる。
!
!javascriptによって内側判定をするとこうなる。
!
!var x = 10; //長方形の左上x
!var y = 10; //長方形の左上y
!var w = 20; //長方形のwidth
!var h = 20; //長方形のheight
!var px = 20; //点のx
!var py = 20; //点のy
!
!function gaiseki(ax,ay,bx,by){
! return ax*by-bx*ay;
!}
!
!function pointInCheck(X,Y,W,H,PX,PY){
! var a = gaiseki(-W,0,PX-W-X,PY-Y);
! var b = gaiseki(0,H,PX-X,PY-Y);
! var c = gaiseki(W,0,PX-X,PY-Y-H);
! var d = gaiseki(0,-H,PX-W-X,PY-H-Y);
! if(a<0&&b<0&&c<0&&d<0) return true;
! else return false;
!}
!
!document.write(pointInCheck(x,y,w,h,px,py));
!
!true
!tureだから、内側にある。
Google Sites
Report abuse