!unvc3d10.f90に加えて、複数材料に対応させてみる。
!さらに、斜めになった部材とかにも木材の異方性を設定できるようにした。
!局所座標系の設定をCYLINDRICALからRECTANGULARに変更した。

! 拘束とか載荷の節点群は、座標指定ではなく、Salomeのmeshのグループで指定する。
! なので、Salomeで、Mesh_1を右クリック-->Create Groups from Geometry-->Nodes
! から、geometryで予めcreate groupしたグループを選択して作っておく。
! 材料別の要素群も同様に、meshのグループを作る必要があるが、
! 要素は並び順情報が要るので、Create Groups from Geometry-->Elements で選択する。

! ちなみに、  /#  で 修正必要箇所 に跳べる。

implicit real*8(a-h,o-z)
dimension x(10000000),y(1000000),z(1000000)
dimension e66(6,6),b(3,3),mesh(30000000,10),ngroup(30000000)
dimension t(3,3),alp(4),bt(4),gmm(4)

!---各種入力データ---
!#--入出力ファイル名
open(7,file='ihou.unv')
open(8,file='kai.inp')
!#--荷重
p=6.d-4!荷重(MN)
!#--Mesh_1 の Mesh Information を見て入力
nnode=541620   !2次要素の節点数(Nodes)
nedge=20911    !1D(edges)の数
nface=154144   !2D(faces)の数
nvolume=287205 !3D(volumes)の数
!---sub mesh の Mesh Information を見て入力
!#--各部材の要素数(volumes)と固定軸回りのオイラー角？
nbuzai=4
pi=asin(1.d0)*2.d0
nyouso1=9448 
alp(1)= 0.d0         ; bt(1)= 1.d0/2.d0*pi ; gmm(1)=0.d0
nyouso2=14790
!alp(2)=11.d0/6.d0*pi ; bt(2)= 0.d0         ; gmm(2)=0.d0        
alp(2)=331.763d0/180.d0*pi ; bt(2)= 0.d0         ; gmm(2)=0.d0        
nyouso3=14655
!alp(3)= 1.d0/6.d0*pi ; bt(3)= 0.d0         ; gmm(3)=0.d0        
alp(3)= 28.237d0/180.d0*pi ; bt(3)= 0.d0         ; gmm(3)=0.d0        
nyouso4=22244
alp(4)= 3.d0/2.d0*pi ; bt(4)= 0.d0         ; gmm(4)=0.d0        
nyouso0=226068    !鋼材
!#--境界条件の節点
nkotei1=12     !固定グループの節点数(Nodes)
nkotei2=12     !固定グループの節点数(Nodes)
nsaika=195     !載荷グループの節点数(Nodes)
!#--変位を読み取る節点(あると便利）
nyomi1=8       !読み取りグループの節点数(Nodes)
nyomi2=7       !読み取りグループの節点数(Nodes)
!#--材料定数
!---鋼材か樹脂のヤング率
young=206.d3
!young=2.84d3
!---鋼材か樹脂のポアソン比
poisson=0.3d0
!---木材のヤング率
ezz=14.461d3
exx=ezz/25.0d0
eyy=ezz/25.0d0
!---木材のポアソン比
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
!---せん断弾性係数
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!---

!---inpファイルの節点情報(*NODE)を作る---
write(8,'("*NODE, NSET=Nall")')
do i=1,19
read (7,'()')!ヘッダーの読み飛ばし
end do
!---節点番号と節点座標の読み込みと書き出し
do i=1,nnode
read (7,*) nset,n1,n2,n3
read(7,*) x(i),y(i),z(i)
write(8,'(i10,a,1pe17.9,a,1pe17.9a,1pe17.9)') nset,',',x(i)/1.0d3&
&,',',y(i)/1.0d3,',',z(i)/1.0d3
end do
!---

!---inpファイルの要素情報を読み取る---
read (7,'()') 
read (7,'()') 
read (7,'()') 
!---Edges情報の読み飛ばし
do i=1,nedge
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) n9,n10
end do
!---Faces情報の読み飛ばし
do i=1,nface
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8,n9,n10,n11
end do
!---Volumes情報の読み込みと書き出し、ついでに体積計算
taiseki=0.d0
do i=1,nvolume
read(7,*) ngroup(i),n1,n2,n3,n4,n5
read(7,*) mesh(ngroup(i),1),mesh(ngroup(i),5),mesh(ngroup(i),2),mesh(ngroup(i),6)&
    &,mesh(ngroup(i),3),mesh(ngroup(i),7),mesh(ngroup(i),8),mesh(ngroup(i),9)
read(7,*) mesh(ngroup(i),10),mesh(ngroup(i),4)
!---体積の計算
x1=x(mesh(ngroup(i),1));y1=y(mesh(ngroup(i),1));z1=z(mesh(ngroup(i),1))
x2=x(mesh(ngroup(i),2));y2=y(mesh(ngroup(i),2));z2=z(mesh(ngroup(i),2))
x3=x(mesh(ngroup(i),3));y3=y(mesh(ngroup(i),3));z3=z(mesh(ngroup(i),3))
x4=x(mesh(ngroup(i),4));y4=y(mesh(ngroup(i),4));z4=z(mesh(ngroup(i),4))
taiseki=taiseki+(&
& (x4-x1)*((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)) &
&+(y4-y1)*((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)) &
&+(z4-z1)*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)) &
&)/6
end do
print*,'体積は',taiseki
read (7,'()') 
read (7,'()') 
read (7,'()') 
!---

!---inpファイルの節点群(*NSET)を作る---
!   以下のグループの並び順は、Salomeで sub mesh を作った順に並び替える
!#--載荷するグループの節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=saika")')
do i=1,int(nsaika/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nsaika,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!#--拘束するグループ1の節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=kousoku1")')
do i=1,int(nkotei1/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nkotei1,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!#--拘束するグループ2の節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=kousoku2")')
do i=1,int(nkotei2/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nkotei2,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!   読み取りグループは、あると便利
!#--読み取りグループ1の節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=yomi1")')
do i=1,int(nyomi1/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nyomi1,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!#--読み取りグループ2の節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=yomi2")')
do i=1,int(nyomi2/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nyomi2,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!---

!---inpファイルの要素群(*ELSET)を作る---
nyou=0
!#--材料1
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E1")')
do i=1,int(nyouso1/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso1,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!#--材料2
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E2")')
do i=1,int(nyouso2/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso2,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!#--材料3
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E3")')
do i=1,int(nyouso3/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso3,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!#--材料4
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E4")')
do i=1,int(nyouso4/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso4,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!#--材料0
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E0")')
do i=1,int(nyouso0/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso0,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!---
close(7)

!#--境界条件(拘束方向は、x:1,y:2,z:3)
write(8,'("*BOUNDARY")')
write(8,'("kousoku1,1,3")')
write(8,'("*BOUNDARY")')
write(8,'("kousoku2,1")')
write(8,'("*BOUNDARY")')
write(8,'("kousoku2,3")')

!#--材料の設定0
write(8,'("*MATERIAL,NAME=Ekou")')
write(8,'("*ELASTIC")')
write(8,'(1pe17.9,a,1pe17.9)')young,',',poisson
write(8,'("*SOLID SECTION,ELSET=E0,MATERIAL=Ekou")')
!   ここから、異方性の設定
!#--材料の設定:木材1
nb=1!部材番号
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
write(8,'("*MATERIAL,NAME=Emoku")')
write(8,'("*ELASTIC,TYPE=ORTHO")')
write(8,'(f12.4,a,f12.4,a,f12.4,a,&
&f12.4,a,f12.4,a,f12.4,a,f12.4,&
& a ,f12.4,a)'),&
&e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
&,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
write(8,'(f12.4)'),e66(4,4)
write(8,'("*ORIENTATION,NAME=OR1,SYSTEM=RECTANGULAR")')
write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
write(8,'("*SOLID SECTION,ELSET=E1,MATERIAL=Emoku,ORIENTATION=OR1")')
!#--材料の設定:木材2
nb=2!部材番号
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
write(8,'("*ORIENTATION,NAME=OR2,SYSTEM=RECTANGULAR")')
write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
write(8,'("*SOLID SECTION,ELSET=E2,MATERIAL=Emoku,ORIENTATION=OR2")')
!#--材料の設定:木材3
nb=3!部材番号
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
write(8,'("*ORIENTATION,NAME=OR3,SYSTEM=RECTANGULAR")')
write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
write(8,'("*SOLID SECTION,ELSET=E3,MATERIAL=Emoku,ORIENTATION=OR3")')
!#--材料の設定:木材4
nb=4!部材番号
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
write(8,'("*ORIENTATION,NAME=OR4,SYSTEM=RECTANGULAR")')
write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
write(8,'("*SOLID SECTION,ELSET=E4,MATERIAL=Emoku,ORIENTATION=OR4")')

!---計算の設定(*STEP~*END STEP)
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
!#--荷重条件(載荷方向は、x:1,y:2,z:3)
write(8,'("*CLOAD")')
write(8,'(A,a,i6,a,1pe25.13)')'saika',',',3,',',-p/nsaika
!#--出力の設定
write(8,'("*NODE PRINT,NSET=yomi1")')
write(8,'("U")')
write(8,'("*NODE PRINT,NSET=yomi2")')
write(8,'("U")')
!write(8,'("*EL PRINT,GLOBAL=YES,ELSET=Eall")')
!write(8,'("S")')
write(8,'("*NODE FILE")')
write(8,'("U")')
write(8,'("*EL FILE,GLOBAL=YES")')
write(8,'("S")')
write(8,'("*END STEP")')
!---
close(8)
end

!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3),t(3,3),alp(4),bt(4),gmm(4)
!応力－ひずみ行列：e66
do i=1,6
do j=1,6
e66(i,j)=0.d0
end do
end do
  
e66(1,1)=1.d0/exx
e66(1,2)=-poixy/exx
e66(1,3)=-poixz/exx
e66(2,1)=-poiyx/eyy
e66(2,2)=1.d0/eyy
e66(2,3)=-poiyz/eyy
e66(3,1)=-poizx/ezz
e66(3,2)=-poizy/ezz
e66(3,3)=1.d0/ezz

!****** 逆行列を取って、応力－ひずみ行列に変換する******
det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
&+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
&-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
do i=1,3
do j=1,3
e66(i,j)=b(j,i)/det
end do
end do
!****** せん断補正係数を入力 ******
!print*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!print*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!print*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
!---固定軸回りのオイラー角？の回転行列
t=0.d0
t(1,1)=              cos(bt(nb))*cos(gmm(nb))
t(2,1)= sin(alp(nb))*sin(bt(nb))*cos(gmm(nb))+cos(alp(nb))*sin(gmm(nb))
t(3,1)=-cos(alp(nb))*sin(bt(nb))*cos(gmm(nb))+sin(alp(nb))*sin(gmm(nb))
t(1,2)=-cos(bt(nb))*sin(gmm(nb))                               
t(2,2)=-sin(alp(nb))*sin(bt(nb))*sin(gmm(nb))+cos(alp(nb))*cos(gmm(nb))
t(3,2)= cos(alp(nb))*sin(bt(nb))*sin(gmm(nb))+sin(alp(nb))*cos(gmm(nb))
t(1,3)= sin(bt(nb))                                                
t(2,3)=-sin(alp(nb))*cos(bt(nb))                                    
t(3,3)= cos(alp(nb))*cos(bt(nb))

return
end
