Commit 6d146bb2 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Create inclusion folder structure

parent 3b83da86

build/inclusion/.gitkeep

0 → 100644
+0 −0

Empty file added.

src/inclusion/.gitkeep

0 → 100644
+0 −0

Empty file added.

src/inclusion/edfb_inclu.f

0 → 100644
+336 −0
Original line number Diff line number Diff line 
      PROGRAM EDFB

CCC   160630

CCC   IES=1 FOR SURROUNDING EXTERNAL SPHERE CENTERED AT ORIGIN;

CCC   HOMOGENEOUS MATERIAL WITHIN THE EXT. SPHERE IS READ IN 

CCC   AS THAT OF THE (NSHL+1)-TH LAYER OF THE 1-ST SPHERE:

CCC   RCF(1,NSHL)=1.0D0, RCF(1,NSHL+1)>1.0D0 AND

CCC   ROS(1) IS RADIUS OF SPHERE 1

CCC

CCC   READ DATA FOR BUILDING VECTOR XIV FROM WITHIN SUB INXI

CCC   

CCC   IDFC>0 WHEN ALL DIEL. FUNCT. OF SPHERES ARE CONSTANTS;

CCC   IDFC=0 WHEN DIEL. FUNCT. OF SPHERES DEPEND ON XI;

CCC   IDFC<0 WHEN DIEL. FUNCT. OF SPHERES ARE AT XIP VALUE ONLY AND

CCC   XI IS SCALE FACTOR FOR DIMENSIONS

CCC   INSN CHOOSES THE VARIABLE THE DIEL. FUNCT. DEPEND ON:

CCC   (INSN=1)=XI;

CCC   (INSN=2)=WN (THE WAVENUMBER, IN m**-1);

CCC   (INSN=3)=WL (THE WAVELENGTH, IN m);

CCC   (INSN=4)=PU (THE ANGULAR FREQUENCY, IN s**-1);

CCC   (INSN=5)=EV (THE PHOTON ENERGY, IN ev);

CCC   INSTPC>0 WHEN VARIABLE INCREASES WITH A CONSTANT STEP;

CCC   INSTPC=0 WHEN VARIABLE IS SAVED IN A VECTOR YOU READ BEFORE

CCC   SUPPORTS EXPERIMENTAL DIELECTRIC FUNCTIONS ONLY

CCC   NSPH=6

      IMPLICIT REAL*8(A-H,O-Z)

CCC   COMMON/C1/DC0(NSHL-NTL+1),DC0M(NSHL-NTL+1,NSPH,NXI),

CCC  1ROS(NSPH),RCF(NSPH,NSHL+1),IOG(NSPH),NSHL(NSPH)

      COMMON/C1/DC0(5),DC0M(5,6,200),

     1ROS(6),RCF(6,9),IOG(6),NSHL(6)

      COMPLEX*16 DC0,DC0M

CCC   COMMON/C3/XIV(NXI),WNS(NXI),WLS(NXI),PUS(NXI),EVS(NXI),

CCC  1VSS(NXI),VNS(5)

      COMMON/C3/XIV(200),WNS(200),WLS(200),PUS(200),EVS(200),

     1VSS(200),VNS(5)

      CHARACTER*3 VNS

 5010 FORMAT(16I5)

 6005 FORMAT(' SPHERE N.',I4)

 6009 FORMAT(' NONTRANSITION LAYER N.',I2,', SCALE = ',A3)

 6010 FORMAT(I5,1X,1PD12.4,1PD12.4)

      IR=5

      IW=6

      IT=7

CCC

CCC   SIZE(I)=VK*ROS(I) IS SIZE PARAMETER in vacuo

CCC

CCC

CCC   READING OF DIELECTRIC FUNCTIONS DRIVEN BY ICI DEFINED BELOW

CCC

      OPEN(IR,FILE='../../test_data/inclusion/DEDFB',STATUS='OLD')

      READ(IR,*)NSPH,IES

      IF(IES.NE.0)IES=1

      READ(IR,*)EXDC,WP,XIP,IDFC,NXI,INSTPC,INSN

      OPEN(IW,FILE='OEDFB',STATUS='UNKNOWN')

      CALL INXI(IR,IW,WP,XIP,IDFC,NXI,INSTPC,INSN)

      READ(IR,5010)(IOG(I),I=1,NSPH)

      DO 113 I=1,NSPH

      IF(IOG(I).LT.I)GO TO 113

      READ(IR,*)NSHL(I),ROS(I)

      NSH=NSHL(I)

      IF(I.EQ.1)NSH=NSH+IES

      DO 112 NS=1,NSH

  112 READ(IR,*)RCF(I,NS)

  113 CONTINUE

      OPEN(IT,FILE='TEDF',FORM='UNFORMATTED',STATUS='UNKNOWN')

      WRITE(IT)NSPH

      WRITE(IT)(IOG(I),I=1,NSPH)

      WRITE(IT)EXDC,WP,XIP,IDFC,NXI

      WRITE(IT)(XIV(I),I=1,NXI)

      DO 115 I=1,NSPH

      IF(IOG(I).LT.I)GO TO 115

      WRITE(IT)NSHL(I),ROS(I)

      NSH=NSHL(I)

      IF(I.EQ.1)NSH=NSH+IES

      WRITE(IT)(RCF(I,NS),NS=1,NSH)

  115 CONTINUE

      DO 468 JXI=1,NXI

      IF((IDFC.NE.0).AND.(JXI.GT.1))GO TO 468

      DO 162 I=1,NSPH

      IF(IOG(I).LT.I)GO TO 162

CCC

      NSH=NSHL(I)

      ICI=(NSH+1)/2

      IF(I.EQ.1)ICI=ICI+IES

      DO 157 IC=1,ICI

      READ(IR,*)DC0(IC)

  157 DC0M(IC,I,JXI)=DC0(IC)

CCC

      WRITE(IT)(DC0(IC),IC=1,ICI)

  162 CONTINUE

  468 CONTINUE

      IF(IDFC.EQ.0)GO TO 474

      WRITE(IW,*)' DIELECTRIC CONSTANTS'

      DO 473 I=1,NSPH

      IF(IOG(I).NE.I)GO TO 473

      ICI=(NSHL(I)+1)/2

      IF(I.EQ.1)ICI=ICI+IES

      WRITE(IW,6005)I

      DO 472 IC=1,ICI

      WRITE(IW,6010)IC,DC0M(IC,I,1)

  472 CONTINUE

  473 CONTINUE

      GO TO 499

  474 WRITE(IW,*)' DIELECTRIC FUNCTIONS'

      DO 478 I=1,NSPH

      IF(IOG(I).NE.I)GO TO 478

      ICI=(NSHL(I)+1)/2

      IF(I.EQ.1)ICI=ICI+IES

      WRITE(IW,6005)I

      DO 477 IC=1,ICI

      WRITE(IW,6009)IC,VNS(INSN)

      DO 476 JXI=1,NXI

      WRITE(IW,6010)JXI,DC0M(IC,I,JXI)

  476 CONTINUE

  477 CONTINUE

  478 CONTINUE

  499 CLOSE(IR)

      CLOSE(IW)

      CLOSE(IT)

      STOP

      END

      SUBROUTINE INXI(IR,IW,WP,XIP,IDFC,NXI,INSTPC,INSN)

      IMPLICIT REAL*8(A-H,O-Z)

      COMMON/C3/XIV(200),WNS(200),WLS(200),PUS(200),EVS(200),

     1VSS(200),VNS(5)

      CHARACTER*3 VNS

      PIGT=DACOS(0.0D0)*4.0D0

      EVC=6.5821188D-16

      IF(IDFC.LT.0)GO TO 300

      IF(INSTPC.EQ.0)GO TO 200

CCC   VLST=V+(NXI-1)*VSTP

      GO TO(105,125,145,165,185),INSN

      RETURN

  105 READ(IR,*)XI,XISTP

      DO 110 JXI=1,NXI

      PU=XI*WP

      WN=PU/3.0D08

      VNS(INSN)='XIV'

      VSS(JXI)=XI

      XIV(JXI)=XI

      PUS(JXI)=PU

      EVS(JXI)=PU*EVC

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

  110 XI=XI+XISTP

      WRITE(IW,6601)

 6601 FORMAT

     1(2X,'JXI',5X,'XIV',10X,'WNS',10X,'WLS',10X,'PUS',10X,'EVS')

      WRITE(IW,6600)

     1(JXI,XIV(JXI),WNS(JXI),WLS(JXI),PUS(JXI),EVS(JXI),JXI=1,NXI)

 6600 FORMAT((I5,5(1PD13.4)))

      RETURN

  125 READ(IR,*)WN,WNSTP

      DO 130 JXI=1,NXI

      XI=3.0D08*WN/WP

      PU=XI*WP

      VNS(INSN)='WNS'

      VSS(JXI)=WN

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

      XIV(JXI)=XI

      PUS(JXI)=PU

      EVS(JXI)=PU*EVC

  130 WN=WN+WNSTP

      WRITE(IW,6602)

 6602 FORMAT

     1(2X,'JXI',5X,'WNS',10X,'WLS',10X,'PUS',10X,'EVS',10X,'XIV')

      WRITE(IW,6600)

     1(JXI,WNS(JXI),WLS(JXI),PUS(JXI),EVS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  145 READ(IR,*)WL,WLSTP

      DO 150 JXI=1,NXI

      WN=PIGT/WL

      XI=3.0D08*WN/WP

      PU=XI*WP

      VNS(INSN)='WLS'

      VSS(JXI)=WL

      WLS(JXI)=WL

      WNS(JXI)=WN

      XIV(JXI)=XI

      PUS(JXI)=PU

      EVS(JXI)=PU*EVC

  150 WL=WL+WLSTP

      WRITE(IW,6603)

 6603 FORMAT

     1(2X,'JXI',5X,'WLS',10X,'WNS',10X,'PUS',10X,'EVS',10X,'XIV')

      WRITE(IW,6600)

     1(JXI,WLS(JXI),WNS(JXI),PUS(JXI),EVS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  165 READ(IR,*)PU,PUSTP

      DO 170 JXI=1,NXI

      XI=PU/WP

      WN=PU/3.0D08

      VNS(INSN)='PUS'

      VSS(JXI)=PU

      PUS(JXI)=PU

      XIV(JXI)=XI

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

      EVS(JXI)=PU*EVC

  170 PU=PU+PUSTP

      WRITE(IW,6604)

 6604 FORMAT

     1(2X,'JXI',5X,'PUS',10X,'WNS',10X,'WLS',10X,'EVS',10X,'XIV')

      WRITE(IW,6600)

     1(JXI,PUS(JXI),WNS(JXI),WLS(JXI),EVS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  185 READ(IR,*)EV,EVSTP

      DO 190 JXI=1,NXI

      PU=EV/EVC

      XI=PU/WP

      WN=PU/3.0D08

      VNS(INSN)='EVS'

      VSS(JXI)=EV

      EVS(JXI)=EV

      PUS(JXI)=PU

      XIV(JXI)=XI

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

  190 EV=EV+EVSTP

      WRITE(IW,6605)

 6605 FORMAT

     1(2X,'JXI',5X,'EVS',10X,'WNS',10X,'WLS',10X,'PUS',10X,'XIV')

      WRITE(IW,6600)

     1(JXI,EVS(JXI),WNS(JXI),WLS(JXI),PUS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  200 GO TO(205,225,245,265,285),INSN

      RETURN

  205 DO 210 JXI=1,NXI

      READ(IR,*)XI

      PU=XI*WP

      WN=PU/3.0D08

      VNS(INSN)='XIV'

      VSS(JXI)=XI

      XIV(JXI)=XI

      PUS(JXI)=PU

      EVS(JXI)=PU*EVC

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

  210 CONTINUE

      WRITE(IW,6601)

      WRITE(IW,6600)

     1(JXI,XIV(JXI),WNS(JXI),WLS(JXI),PUS(JXI),EVS(JXI),JXI=1,NXI)

      RETURN

  225 DO 230 JXI=1,NXI

      READ(IR,*)WN

      XI=3.0D08*WN/WP

      PU=XI*WP

      VNS(INSN)='WNS'

      VSS(JXI)=WN

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

      XIV(JXI)=XI

      PUS(JXI)=PU

      EVS(JXI)=PU*EVC

  230 CONTINUE

      WRITE(IW,6602)

      WRITE(IW,6600)

     1(JXI,WNS(JXI),WLS(JXI),PUS(JXI),EVS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  245 DO 250 JXI=1,NXI

      READ(IR,*)WL

      WN=PIGT/WL

      XI=3.0D08*WN/WP

      PU=XI*WP

      VNS(INSN)='WLS'

      VSS(JXI)=WL

      WLS(JXI)=WL

      WNS(JXI)=WN

      XIV(JXI)=XI

      PUS(JXI)=PU

      EVS(JXI)=PU*EVC

  250 CONTINUE

      WRITE(IW,6603)

      WRITE(IW,6600)

     1(JXI,WLS(JXI),WNS(JXI),PUS(JXI),EVS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  265 DO 270 JXI=1,NXI

      READ(IR,*)PU

      XI=PU/WP

      WN=PU/3.0D08

      VNS(INSN)='PUS'

      VSS(JXI)=PU

      PUS(JXI)=PU

      XIV(JXI)=XI

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

      EVS(JXI)=PU*EVC

  270 CONTINUE

      WRITE(IW,6604)

      WRITE(IW,6600)

     1(JXI,PUS(JXI),WNS(JXI),WLS(JXI),EVS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  285 DO 290 JXI=1,NXI

      READ(IR,*)EV

      PU=EV/EVC

      XI=PU/WP

      WN=PU/3.0D08

      VNS(INSN)='EVS'

      VSS(JXI)=EV

      EVS(JXI)=EV

      PUS(JXI)=PU

      XIV(JXI)=XI

      WNS(JXI)=WN

      WLS(JXI)=PIGT/WN

  290 CONTINUE

      WRITE(IW,6605)

      WRITE(IW,6600)

     1(JXI,EVS(JXI),WNS(JXI),WLS(JXI),PUS(JXI),XIV(JXI),JXI=1,NXI)

      RETURN

  300 IF(INSTPC.GT.0)GO TO 315

      DO 310 JXI=1,NXI

      READ(IR,*)XI

      VNS(INSN)='XIV'

      VSS(JXI)=XI

      XIV(JXI)=XI

  310 CONTINUE

      GO TO 330

  315 READ(IR,*)XI,XISTP

      DO 320 JXI=1,NXI

      VNS(INSN)='XIV'

      VSS(JXI)=XI

      XIV(JXI)=XI

  320 XI=XI+XISTP

  330 PU=XIP*WP

      WN=PU/3.0D08

      WRITE(IW,6611)

 6611 FORMAT

     1(10X,'XIP',10X,'WN ',10X,'WL ',10X,'PU ',10X,'EV')

      WRITE(IW,6610)XIP,WN,PIGT/WN,PU,PU*EVC

 6610 FORMAT((5X,5(1PD13.4)))

      WRITE(IW,*)' SCALE FACTORS XI'

      WRITE(IW,6612)(JXI,XIV(JXI),JXI=1,NXI)

 6612 FORMAT(I5,1PD13.4)

      RETURN

      END

CCC

src/inclusion/inclu.f

0 → 100644
+3335 −0

File added.

Preview size limit exceeded, changes collapsed.

test_data/inclusion/.gitkeep

0 → 100644
+0 −0

Empty file added.