IMSL를 이용한 BSPLINE

---

공학문제를 풀다보면 보간이 필요 한 경우가 있는데, 

1) 이 때 IMSL를 이용하여 이를 계산하고자 한다. 

2) 그리고 이 후 이를 C#을 이용하여 GUI환경의 프로그램으로 만들고자 한다.

3) 그 결과 일반적인 airfoil data를 수정한다.

4) 또한 airfoil data를 ICEM-CFD로 읽고, LES용 mesh를 생성한다. 

5) 그리고 광대역 뒷전 소음을 예측을 위해 Patran에서 airfoil surface mesh 생성한다.

오늘은 그 첫 시간으로 IMSL을 이용한 DATA 보간을 수행하며, 이 때 BSpline 방법이 이용된다.




!============================================
!                    CODE

!============================================

!DEC$ OBJCOMMENT LIB:'LIBGUIDE.LIB'

PROGRAM BSPLINE_INTERPOLATION

INCLUDE 'LINK_FNL_STATIC.H'

USE SPLINE_FITTING_INT

USE SHOW_INT

USE NORM_INT

IMPLICIT NONE

! THIS IS EXAMPLE 1 FOR SPLINE_FITTING, NATURAL SPLINE

! INTERPOLATION USING CUBIC SPLINES. USE THE FUNCTION

! EXP(-X**2/2) TO GENERATE SAMPLES.

INTEGER :: I, NI

INTEGER, PARAMETER :: NDATA=24, NORD=4, NDEGREE=NORD-1,NBKPT=NDATA+2*NDEGREE, NCOEFF=NBKPT-NORD, NVALUES=3*NDATA


! NDATA: 보간 할 원본 데이터의 개수
! NVALUES : 보간 된 출력 데이터의 개수

REAL(KIND(1E0)), PARAMETER :: ZERO=0E0, ONE=1E0, HALF=5E-1

REAL(KIND(1E0)), PARAMETER :: DELTA_X=0.15, DELTA_XV=0.4*DELTA_X

REAL(KIND(1E0)), TARGET :: XDATA(NDATA), YDATA(NDATA), SPLINE_DATA (3, NDATA), BKPT(NBKPT)

REAL(KIND(1E0)), TARGET :: YCHECK(NVALUES), COEFF(NCOEFF), XVALUES(NVALUES), YVALUES(NVALUES), DIFFS

REAL(KIND(1E0)), POINTER :: POINTER_BKPT(:)

TYPE (S_SPLINE_KNOTS) BREAK_POINTS

TYPE (S_SPLINE_CONSTRAINTS) CONSTRAINTS(2)

!============INPUT VALUE===========
! 원본 데이터를 계산한다.

XDATA = (/((I-1)*DELTA_X, I=1,NDATA)/)          !ORIGINAL XDATA

YDATA = EXP(-HALF*XDATA**2)                     !ORIGINAL YDATA

!=====DISPLAY INPUT VALUE===========

DO NI=1, NDATA

    WRITE(*,*) XDATA(NI), YDATA(NI)   

ENDDO

XVALUES =(/(0.03+(I-1)*DELTA_XV,I=1,NVALUES)/)      !SHIFTED ORIGINAL XDATA

YCHECK= EXP(-HALF*XVALUES**2)                       !SHIFTED ORIGINAL YDATA

SPLINE_DATA(1,:)=XDATA                              !IN/OUTPUT VARIABLE

SPLINE_DATA(2,:)=YDATA                              !IN/OUTPUT VARIABLE

SPLINE_DATA(3,:)=ONE

! DEFINE THE KNOTS FOR THE INTERPOLATION PROBLEM.

BKPT(1:NDEGREE) = (/(I*DELTA_X, I=-NDEGREE,-1)/)

BKPT(NORD:NBKPT-NDEGREE) = XDATA

BKPT(NBKPT-NDEGREE+1:NBKPT) = (/(XDATA(NDATA)+I*DELTA_X, I=1,NDEGREE)/)

! ASSIGN THE DEGREE OF THE POLYNOMIAL AND THE KNOTS.

POINTER_BKPT => BKPT

BREAK_POINTS=S_SPLINE_KNOTS(NDEGREE, POINTER_BKPT)

! THESE ARE THE NATURAL CONDITIONS FOR INTERPOLATING CUBIC

! SPLINES. THE DERIVATIVES MATCH THOSE OF THE INTERPOLATING

! FUNCTION AT THE ENDS.

CONSTRAINTS(1)=SPLINE_CONSTRAINTS (DERIVATIVE=2, POINT=BKPT(NORD), TYPE='==', VALUE=-ONE)

CONSTRAINTS(2)=SPLINE_CONSTRAINTS (DERIVATIVE=2,POINT=BKPT(NBKPT-NDEGREE), TYPE= '==',VALUE=(-ONE+XDATA(NDATA)**2)*YDATA(NDATA))

COEFF = SPLINE_FITTING(DATA=SPLINE_DATA, KNOTS=BREAK_POINTS,CONSTRAINTS=CONSTRAINTS)

YVALUES=SPLINE_VALUES(0, XVALUES, BREAK_POINTS, COEFF)

DIFFS=NORM(YVALUES-YCHECK,HUGE(1))/DELTA_X**NORD

!=====DISPLAY OUTPUT VALUE===========

DO NI=1, NVALUES

    IF ( XVALUES(NI) <= XDATA(NDATA))        WRITE(*,*) XVALUES(NI), YVALUES(NI)   

ENDDO

IF (DIFFS <= ONE) THEN

    WRITE(*,*) 'EXAMPLE 1 FOR SPLINE_FITTING IS CORRECT.'

END IF

STOP

ENDPROGRAM BSPLINE_INTERPOLATION

!============================================

!                   Result

!============================================

  0.0000000E+00   1.000000

  0.1500000      0.9888130

  0.3000000      0.9559975

  0.4500000      0.9037071

  0.6000000      0.8352702

  0.7500000      0.7548396

  0.9000000      0.6669768

   1.050000      0.5762290

   1.200000      0.4867522

   1.350000      0.4020213

   1.500000      0.3246525

   1.650000      0.2563401

   1.800000      0.1978987

   1.950000      0.1493818

   2.100000      0.1102505

   2.250000      7.9559512E-02

   2.400000      5.6134757E-02

   2.550000      3.8725752E-02

   2.700000      2.6121404E-02

   2.850000      1.7227467E-02

   3.000000      1.1108996E-02

   3.150000      7.0041651E-03

   3.300000      4.3178373E-03

   3.450000      2.6025851E-03

  2.9999999E-02  0.9995434

  9.0000004E-02  0.9959494

  0.1500000      0.9888129

  0.2100000      0.9781896

  0.2700000      0.9642053

  0.3300000      0.9470045

  0.3900000      0.9267669

  0.4500000      0.9037071

  0.5100000      0.8780502

  0.5700000      0.8500577

  0.6300000      0.8200007

  0.6900000      0.7881626

  0.7500000      0.7548395

  0.8100000      0.7203276

  0.8700000      0.6849223

  0.9300000      0.6489181

  0.9900000      0.6125971

   1.050000      0.5762289

   1.110000      0.5400755

   1.170000      0.5043683

   1.230000      0.4693307

   1.290000      0.4351585

   1.350000      0.4020213

   1.410000      0.3700771

   1.470000      0.3394438

   1.530000      0.3102281

   1.590000      0.2825097

   1.650000      0.2563401

   1.710000      0.2317623

   1.770000      0.2087859

   1.830000      0.1874117

   1.890000      0.1676222

   1.950000      0.1493817

   2.010000      0.1326495

   2.070000      0.1173673

   2.130000      0.1034725

   2.190000      9.0895399E-02

   2.250000      7.9559527E-02

   2.310000      6.9387339E-02

   2.370000      6.0298394E-02

   2.430000      5.2211456E-02

   2.490000      4.5046724E-02

   2.550000      3.8725760E-02

   2.610000      3.3171590E-02

   2.670000      2.8312404E-02

   2.730000      2.4077961E-02

   2.790000      2.0403031E-02

   2.850000      1.7227462E-02

   2.910000      1.4493150E-02

   2.970000      1.2149527E-02

   3.030000      1.0148127E-02

   3.090000      8.4457034E-03

   3.150000      7.0041865E-03

   3.210000      5.7872944E-03

   3.270000      4.7649527E-03

   3.330000      3.9088321E-03

   3.390000      3.1947396E-03

   3.450000      2.6026356E-03

Press any key to continue . . .

---