* PACKAGE ASETS  !" ���� ���(���ɸ��)
*
*"  [HIS] 92/09/17(takepiro)
*"        93/04/30(takepiro) ���ե�����̾���ѹ�
*"        93/05/27(takepiro) �Ť߷׻�
*
***********************************************************************
      SUBROUTINE SETRAD      !"  ���٥� 
     O         ( ARAD  , DRAD  , HARAD ,
     D           KDIM                    )
*
*"               (* �����ӥ����եݥ���� ���ɸ *)
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCCOM)                      !" ɸ��ʪ�����
      INCLUDE   (ZHDIM)                      !" ʸ����ʸ����
#else
#include        "zccom.F"                    !" ɸ��ʪ�����
#include        "zhdim.F"                    !" ʸ����ʸ����
#endif
      INTEGER    KDIM
*
*   [OUTPUT] 
      REAL       ARAD  ( 0:KDIM )              !" ���٥�
      REAL       DRAD  ( 0:KDIM )              !" ����
      CHARACTER  HARAD  *(*)                 !" ��ľ��̾��
*
*    [INTERNAL WORK]    
      INTEGER    KDMAX
      PARAMETER  ( KDMAX = 100 )
      REAL       XSI  ( 0:KDMAX )              !" �Υ�٥�
      REAL       DXSI ( 0:KDMAX )              !" ����
      INTEGER    K
*
      WRITE ( 6,* ) ' RADIAL LEVEL RAD A0.4 DATE=93/05/27'
*
*"         < 1. ̾�Υ��å� >
*
      CALL STNRAD
     O         ( HARAD ,
     I           KDIM   )
*
*"         < 2. ��, ����η׻� >
*
      CALL SETXSI
     O         ( XSI , DXSI , 
     I           KDIM         )
*
      DO 3100 K = 0 , KDIM
        ARAD  ( K ) = ( ( ERI + ERO ) + XSI(K) * ( ERO - ERI ) )*0.5
 3100 CONTINUE
*
      DO 3200 K = 0, KDIM
        DRAD  ( K ) = DXSI( K ) * ( ERO - ERI ) * 0.5 * ARAD(K)**2
 3200 CONTINUE
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SETXSI     !" �����ӥ����եݥ���Ȥη׻�
     O              (  XSI , DXSI , 
     I                 KDIM         )
*
*   [PARAM] 
      INTEGER    KDIM
*
*   [OUTPUT] 
      REAL       XSI  ( 0:KDIM )               !" �Υ�٥�
      REAL       DXSI ( 0:KDIM )               !" �Υ�٥�
*
*    [INTERNAL WORK]    
      REAL       PI                          !" �߼�Ψ
      INTEGER    J, K
*
*"  < 1. �����ӥ����ե�å���ݥ���Ȥη׻� > 
*
      PI = ATAN( 1.0D0 ) * 4.0D0
*
      DO 1000 K = 0, KDIM
         XSI ( K ) = COS ( PI * K / KDIM )
 1000 CONTINUE 
*
*"  < 2. �Ťߤη׻� > 
*
      DO 2100 K = 0, KDIM
         DXSI( K ) = 0. 
         DO 2200 J = 0, KDIM, 2
            DXSI( K ) = DXSI( K ) 
     &                  + 2. / (1-J**2) * COS ( J*K*PI / KDIM )
 2200    CONTINUE 
         DXSI( K ) = DXSI( K ) 
     &               - 1 
     &              -  0.5 * ( (-1)**(KDIM-1) -1 ) / ( 1 - KDIM**2 )
         DXSI( K ) = 2. / KDIM * DXSI( K )
 2100 CONTINUE 
*
      DXSI( 0 )    =  DXSI( 0 ) / 2 
      DXSI( KDIM ) =  DXSI( KDIM ) / 2 
*
      RETURN
      END
***********************************************************************
      SUBROUTINE STNRAD
     O               ( HARAD ,
     I                 KDIM   )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZHDIM)                      !" ʸ����ʸ����
      INCLUDE   (ZCCOM)                      !" ɸ��ʪ�����
#else
#include        "zhdim.F"                    !" ʸ����ʸ����
#include        "zccom.F"                    !" ɸ��ʪ�����
#endif
*
      INTEGER    KDIM
*
*   [OUTPUT] 
      CHARACTER  HARAD  *(*)                 !" ��ľ��̾�Ρ�����
*
*    [INTERNAL WORK]    
      INTEGER    INUM, IAZ
      CHARACTER  HNUM   *(NCC)
*
      CHARACTER  HARAD1 *(NCC)
      DATA       HARAD1 / 'RAD'  /
*
*    [EXTERNAL FUNC]
      INTEGER    LENC
*
*"         < 1. ̾�Υ��å� >
*
      CALL GULCHR
     I     ( '(I4)', KDIM  ,
     O       HNUM  , INUM   )
*
      IAZ    = LENC( HARAD1 )
      HARAD  = HARAD1(1:IAZ)//'P'//HNUM(1:INUM)
*
      WRITE ( HNUM, '(F6.2)' ) ETA
      CALL CLADJ ( HNUM )
      INUM = LENC( HNUM )
*
      IAZ    = LENC( HARAD )
      HARAD  = HARAD(1:IAZ)//'-'//HNUM(1:INUM)
*
      RETURN
      END
