Programmer 7500

home *** CD-ROM | disk | FTP | other *** search

/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / PROGJRN / PJ_6_6.ARC / FCODE.ARC / SLITE3.FOR < prev next >

Wrap

Text File | 1987-09-07 | 75.6 KB | 2,416 lines

PROGRAM LITE REAL X0(4,30,7),BETA(4,30,2),PSI(4,30,2),RHO(4,30),XN(30,400,3), 1 DAYF(3,400),DELF(400),F(4,30,30),WSOL(3,400),DSOL(4,30), 2 TAU(10),DW(5),LL(4,30,9),A(4,30),COEF(30,30),LS(3),Y(3), 3 YM(3),XY(3,4,2),LUM(4,30),AK2(10),X(3),COFB(30,1),FS(30), 4 ALPHA(30),FR(30),AIK(30,400),LUMF(30,400),LUMI(30,400), 5 XYW(5,2,2),XSUN(5,4,2),ANG(5),OHW(5),OHREF(5), 6 RING(5),OHTAU(5),RS(4),WW(5),WH(5),WCX(5),WCY(5), 8 W(12),BW(12), 9 PLTTMP(4,3) C INTEGER IDS(30),IOB(4,30,30),JOB(30,30,5),KSUN(4,30),NW(6),TZN, 1 KOB(4,30,30),KD(30,11,11),ICO(30),JCO(30,5),NM(30),NN(30), 2 ISUN(5),NUM(30),IJOB(26),KJOB(4,30,30),KKOB(4,30,30), 3 IOH(5),KOH(5),ISTART(5),ITYPW(5),ISD(5), 4 IMG(5,20,20),ICNR(5,4) INTEGER*4 ISEED C COMMON /BLK1/ X0,LL,XDP,YDP COMMON /BLK2/ NW,PI,IOB,ICLD,ISEED,NBUNDL,A,KJOB COMMON /BLK3/ AK2,TAU,NC COMMON /BLK4/ BETAS,LS,AUX(3),GRNDL COMMON /BLK5/ ICO,JCO,NSW COMMON /BLK6/ XN,AIK,LUMF,LUMI,WSOL,DAYF,DELF,IMG 2 ,ICNR C ISNGRD - DIMENSION OF IMG(), DIRECT SUN IMAGE ARRAY ISNGRD=20 PI=ATAN(1.)*4. RAD=PI/180. NBUNDL=2000 ISEED=123456 C INID=1.. SIMPLIFIED INPUT FOR RECTANGULAR ROOMS C IWRITE=1.. ILLUMINATION DATA FOR WORKING SURFACE ARE PRINTED C =2.. DAYLIGHT FACTORS ARE PRINTED C =3.. ILLUMINATION AND DAYLIGHT FACTORS ARE PRINTED C INID=2.. FULL INPUT FOR ARBITRARY ROOMS, C IWRITE=1.. BRIEF OUTPUT (DATA FOR WORKING SURFACE ONLY) C =2.. DETAILED OUTPUT (DATA FOR ALL SURFACES) C IPLOT=0.. NO PLOT IS BEING GENERATED C =1.. ILLUMINATION LEVELS ARE PLOTTED (EQUALLY-SPACED LINES) C =2.. ILLUMINATION LEVELS ARE PLOTTED (EQUALLY-SPACED LEVELS) C =3.. DAYLIGHT FACTORS ARE PLOTTED (EQUALLY-SPACED LEVELS) C NUMIT= # OF ITERATIONS C IOH=0.. NO OVERHANGS, C =1.. OVERHANG ON TOP, C =2.. OVERHANGS ON BOTH SIDES, C =3.. OVERHANG ON BOTTOM ONLY, C =4.. OVERHANGS ON TOP AND BOTH SIDES, C =5.. OVERHANGS ON ALL FOUR SIDES. C IWMOD=1.. INPUT IS SUN POSITION, SOLAR RADIATION, ETC., C =2.. INPUT IS LATITUDE,LONGITUDE,TIME,ETC. C ITYPW=1.. CLEAR WINDOW, C =2.. WINDOW WITH SHADING DEVICE, C =3.. DIFFUSING WINDOW. C ISD =1.. SHEER CURTAIN C C INPUT ROUTINES OPEN(UNIT=5,FILE='SLITE.IN',STATUS='OLD') OPEN(UNIT=7,FILE='SLITE.OUT',STATUS='NEW') C * RECL=200) OPEN(UNIT=8,FILE='SLITE.PLT',STATUS='NEW') C CALL LOGO(6) READ(5,*) INID,IWMOD,NUMIT,IWRITE,IPLOT IF(INID.NE.1) GOTO 31 C C *** SIMPLIFIED INPUT - LIMITED TO RECTANGULAR ROOM WITH C *** RECTANGULAR WINDOWS, AND SIMPLE OUTSIDE OBSTRUCTIONS C C ROOM SIZE AND LOCATION C READ(5,*) RW,RD,RH READ(5,*) RE,RO C C WINDOW DATA C NWINDO = NUMBER OF WINDOWS, C LOCW=1.. WINDOWS ARE IN THE FRONT WALL, C =2.. WINDOWS ARE IN THE CEILING (SKYLIGHTS) READ(5,*) NWINDO,LOCW C INITIALIZATION OF AUXILIARY VALUES IF (NWINDO.GE.1.AND.NWINDO.LE.5) GOTO 4010 WRITE(7,8002) STOP 4010 NSS=0 NC=0 NSW=0 ND=0 DO 4000 I=1,NWINDO IDS(I)=0 OHW(I)=0. RING(I)=0. OHREF(I)=0. KOH(I)=0 4000 OHTAU(I)=0. C DATA FOR EACH WINDOW DO 4100 I=1,NWINDO READ(5,*) WW(I),WH(I),WCX(I),WCY(I),NM(I),NN(I),RHO(1,I) READ(5,*) ITYPW(I),AK2(I),TAU(I),DW(I),IOH(I) NC=NC+1/ITYPW(I) NSW=NSW+ITYPW(I)/2-ITYPW(I)/3 ND=ND+ITYPW(I)/3 IF (LOCW.EQ.1) IDS(I)=2 ALPHA(I)=1. ISD(I)=1 IF(ITYPW(I).NE.2) GOTO 4050 C SHEER CURTAIN READ(5,*) ALPHA(I) 4050 IF(ITYPW(I).EQ.3) ALPHA(I)=0. IF(IOH(I).EQ.0) GOTO 4100 C OVERHANG DATA READ(5,*) OHW(I),RING(I),OHREF(I),KOH(I),OHTAU(I) IF (KOH(I).EQ.0) OHTAU(I)=0. 4100 CONTINUE C WINDOW ORDERING (PUT IN SEQUENCE CLEAR--SHADED--DIFFUSING) IF(ITYPW(NWINDO).GE.ITYPW(1).OR.NWINDO.GT.2) GOTO 4110 AAA=AK2(1) AK2(1)=AK2(2) AK2(2)=AAA AAA=TAU(1) TAU(1)=TAU(2) TAU(2)=AAA AAA=DW(1) DW(1)=DW(2) DW(2)=AAA AAA=ALPHA(1) ALPHA(1)=ALPHA(2) ALPHA(2)=AAA AAA=OHW(1) OHW(1)=OHW(2) OHW(2)=AAA AAA=RING(1) RING(1)=RING(2) RING(2)=AAA AAA=OHREF(1) OHREF(1)=OHREF(2) OHREF(2)=AAA AAA=OHTAU(1) OHTAU(1)=OHTAU(2) OHTAU(2)=AAA AAA=WW(1) WW(1)=WW(2) WW(2)=AAA AAA=WH(1) WH(1)=WH(2) WH(2)=AAA AAA=WCX(1) WCX(1)=WCX(2) WCX(2)=AAA AAA=WCY(1) WCY(1)=WCY(2) WCY(2)=AAA AAA=RHO(1,1) RHO(1,1)=RHO(1,2) RHO(1,2)=AAA III=ITYPW(1) ITYPW(1)=ITYPW(2) ITYPW(2)=III III=IOH(1) IOH(1)=IOH(2) IOH(2)=III III=KOH(1) KOH(1)=KOH(2) KOH(2)=III III=NM(1) NM(1)=NM(2) NM(2)=III III=NN(1) NN(1)=NN(2) NN(2)=III C C REFLECTIVITIES AND DISCRETIZATION FOR WALLS C C IF LOCW=1, SEQUENCE IS WINDOWS-LEFT WALL-REAR WALL-RIGHT WALL-CEILING- C FLOOR-FRONT WALL; IF LOCW=2, FRONT WALL AND CEILING ARE INTERCHANGED C TO PLACE SURFACE WITH CUT-OUT (WINDOWS) LAST 4110 DO 4120 I=1,6 IC=I+NWINDO-1 IF(I.EQ.1) IC=NWINDO+2*(4-LOCW) IF(I.EQ.5) IC=NWINDO+2*(1+LOCW) 4120 READ(5,*) NM(IC),NN(IC),RHO(1,IC) C C WORKING SURFACE C READ(5,*) NM(NWINDO+7),NN(NWINDO+7),WSE C C DESCRIPTION OF OUTSIDES C READ(5,*) IOPP C C OPPOSING BUILDING C IF(IOPP.EQ.1) READ(5,*) HO,WO,XO,DO,RHO(2,3) C C SUBJECT BUILDING C READ(5,*) HS,WS,XS,RHO(2,1) C SET AUXILIARY VALUES NW(2)=IOPP+2 NWC=5 NWO=1 NWS=1 IF (LOCW.EQ.1) NSS=1 N0=NC+NSW N01=NC+1 N1=N0+ND N11=N1+1 N2=N1+NWC NT=N2+NWO NT1=NT+1 NTS=NT1 NSP1=NSS+1 C CALCULATE VALUES FOR CUT-OUT AND OBSTRUCTION IDENTIFIERS IMAX=NWINDO+6 IMAX1=IMAX+1 ICO(IMAX)=NWINDO DO 4170 I=1,IMAX1 DO 4160 J=1,IMAX 4160 IOB(1,I,J)=0 4170 IOB(1,I,I)=-1 DO 4175 I=1,NWINDO JCO(IMAX,I)=I IOB(1,I,IMAX)=-1 4175 IOB(1,IMAX,I)=-1 IOB(1,IMAX1,NWINDO+5)=-1 C DO 4180 I=1,IMAX1 DO 4180 J=1,3 IOB(2,I,J)=0 IOB(2,J,J)=-1 4180 KOB(2,I,J)=-1 DO 4190 I=1,5 DO 4190 J=2,3 4190 KOB(2,I+NWINDO,J)=0 KOB(2,NWINDO+5,2)=-1 KOB(2,NWINDO+7,3)=0 RO=RO*RAD SRO=SIN(RO) CRO=COS(RO) C C DESCRIPTION OF WINDOWS C DO 4200 I=1,NWINDO X0(1,I,1)=(WCX(I)-WW(I)/2.)*SRO X0(1,I,2)=-(WCX(I)-WW(I)/2.)*CRO X0(1,I,3)=RE+WCY(I)-WH(I)/2. X0(1,I,4)=WW(I) X0(1,I,5)=0. X0(1,I,6)=WW(I) X0(1,I,7)=WH(I) BETA(1,I,1)=PI/2. BETA(1,I,2)=0. PSI(1,I,1)=RO-PI/2. PSI(1,I,2)=0. IF (LOCW.NE.2) GOTO 4200 X0(1,I,1)=X0(1,I,1) - (WCY(I)-WH(I)/2.)*CRO X0(1,I,2)=-WCX(I)*CRO-(WCY(I)-WH(I)/2.)*SRO+WW(I)*.5*SRO X0(1,I,3)=RE+RH BETA(1,I,2)=PI/2. PSI(1,I,2)=RO-PI 4200 CONTINUE C C FRONT WALL C I=NWINDO+2*(4-LOCW) X0(1,I,1)=-RW/2.*SRO X0(1,I,2)=RW/2.*CRO X0(1,I,3)=RE X0(1,I,4)=RW X0(1,I,5)=0. X0(1,I,6)=RW X0(1,I,7)=RH BETA(1,I,1)=PI/2. BETA(1,I,2)=0. PSI(1,I,1)=RO-PI/2. PSI(1,I,2)=0. C C LEFT WALL C I=NWINDO+1 X0(1,I,1)=-RD*CRO-RW/2.*SRO X0(1,I,2)=-RD*SRO+RW/2.*CRO X0(1,I,3)=RE X0(1,I,4)=RD X0(1,I,5)=0. X0(1,I,6)=RD X0(1,I,7)=RH BETA(1,I,1)=PI/2. BETA(1,I,2)=0. PSI(1,I,1)=RO PSI(1,I,2)=0. C C REAR WALL C I=NWINDO+2 X0(1,I,1)=-RD*CRO+RW/2.*SRO X0(1,I,2)=-RD*SRO-RW/2.*CRO X0(1,I,3)=RE X0(1,I,4)=RW X0(1,I,5)=0. X0(1,I,6)=RW X0(1,I,7)=RH BETA(1,I,1)=PI/2. BETA(1,I,2)=0. PSI(1,I,1)=RO+PI/2. PSI(1,I,2)=0. C C RIGHT SIDE C I=NWINDO+3 X0(1,I,1)=RW/2.*SRO X0(1,I,2)=-RW/2.*CRO X0(1,I,3)=RE X0(1,I,4)=RD X0(1,I,5)=0. X0(1,I,6)=RD X0(1,I,7)=RH BETA(1,I,1)=PI/2. BETA(1,I,2)=0. PSI(1,I,1)=RO+PI PSI(1,I,2)=0. C C CEILING C I=NWINDO+2*(1+LOCW) X0(1,I,1)=RW/2.*SRO-RD*CRO X0(1,I,2)=-RW/2.*CRO-RD*SRO X0(1,I,3)=RE+RH X0(1,I,4)=RW X0(1,I,5)=0. X0(1,I,6)=RW X0(1,I,7)=RD BETA(1,I,1)=PI/2. BETA(1,I,2)=PI/2. PSI(1,I,1)=RO+PI/2. PSI(1,I,2)=RO C C FLOOR C I=NWINDO+5 X0(1,I,1)=-RW/2.*SRO-RD*CRO X0(1,I,2)=RW/2.*CRO-RD*SRO X0(1,I,3)=RE X0(1,I,4)=RW X0(1,I,5)=0. X0(1,I,6)=RW X0(1,I,7)=RD BETA(1,I,1)=PI/2. BETA(1,I,2)=PI/2. PSI(1,I,1)=RO-PI/2. PSI(1,I,2)=RO C C WORKING SURFACE C I=NWINDO+7 X0(1,I,1)=-RW/2.*SRO-RD*CRO X0(1,I,2)=RW/2.*CRO-RD*SRO X0(1,I,3)=RE+WSE X0(1,I,4)=RW X0(1,I,5)=0. X0(1,I,6)=RW X0(1,I,7)=RD BETA(1,I,1)=PI/2. BETA(1,I,2)=PI/2. PSI(1,I,1)=RO-PI/2. PSI(1,I,2)=RO RHO(1,I)=1. C C OUTSIDE WINDOW WALL C X0(2,1,1)=(.5*WS+XS)*SRO X0(2,1,2)=-(.5*WS+XS)*CRO X0(2,1,3)=0. X0(2,1,4)=WS X0(2,1,5)=0. X0(2,1,6)=WS X0(2,1,7)=HS BETA(2,1,1)=PI/2. BETA(2,1,2)=0. PSI(2,1,1)=RO+PI/2. PSI(2,1,2)=0. C C OUTSIDE GROUND C IF (IOPP.EQ.-1) GOTO 4210 HGR=WS IF(IOPP.EQ.1) HGR=DO X0(2,2,1)=-(0.75*WS-XS)*SRO X0(2,2,2)= (0.75*WS-XS)*CRO X0(2,2,3)=0. X0(2,2,4)=1.5*WS X0(2,2,5)=0. X0(2,2,6)=1.5*WS X0(2,2,7)=HGR BETA(2,2,1)=PI/2. BETA(2,2,2)=PI/2. PSI(2,2,1)=RO-PI/2. PSI(2,2,2)=RO IF(IOPP.LE.0) GOTO 4210 C C OPPOSING BUILDING C X0(2,3,1)=DO*CRO-(.5*WO-XO)*SRO X0(2,3,2)=DO*SRO+(.5*WO-XO)*CRO X0(2,3,3)=0. X0(2,3,4)=WO X0(2,3,5)=0. X0(2,3,6)=WO X0(2,3,7)=HO BETA(2,3,1)=PI/2. BETA(2,3,2)=0. PSI(2,3,1)=RO-PI/2. PSI(2,3,2)=0. C C INPUT DATA PRINT-OUT C 4210 RO=RO/RAD WRITE(7,6300) RW,RD,RH,RO IF(LOCW.EQ.1) WRITE(7,6302) NWINDO IF(LOCW.EQ.2) WRITE(7,6303) NWINDO DO 4220 I=1,NWINDO WRITE(7,6310) I,WW(I),WH(I),WCX(I),WCY(I),NM(I),NN(I) IF(ITYPW(I).EQ.1) WRITE(7,6311) ITYPW(I) IF(ITYPW(I).NE.2) GOTO 4215 IF(ISD(I).EQ.1) WRITE(7,6312) ITYPW(I),ALPHA(I) 4215 IF(ITYPW(I).EQ.3) WRITE(7,6313) ITYPW(I) WRITE(7,6320) AK2(I),DW(I),TAU(I),RHO(1,I) IF(IOH(I).EQ.0) GOTO 4220 IF(IOH(I).EQ.1) WRITE(7,6031) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.2) WRITE(7,6032) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.3) WRITE(7,6033) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.4) WRITE(7,6034) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.5) WRITE(7,6035) OHW(I),RING(I),OHREF(I) IF(KOH(I).EQ.0) WRITE(7,6036) IF(KOH(I).EQ.1) WRITE(7,6037) OHTAU(I) IF(KOH(I).EQ.2) WRITE(7,6038) OHTAU(I) 4220 CONTINUE NWP6=NWINDO+2*(4-LOCW) NWP4=NWINDO+2*(1+LOCW) WRITE(7,6331) NM(NWP6),NN(NWP6),RHO(1,NWP6) WRITE(7,6332) NM(NWINDO+1),NN(NWINDO+1),RHO(1,NWINDO+1) WRITE(7,6333) NM(NWINDO+2),NN(NWINDO+2),RHO(1,NWINDO+2) WRITE(7,6334) NM(NWINDO+3),NN(NWINDO+3),RHO(1,NWINDO+3) WRITE(7,6335) NM(NWP4),NN(NWP4),RHO(1,NWP4) WRITE(7,6336) NM(NWINDO+5),NN(NWINDO+5),RHO(1,NWINDO+5) WRITE(7,6337) NM(NWINDO+7),NN(NWINDO+7),WSE WRITE(7,6340) WS,HS,XS,RHO(2,1) IF(IOPP.EQ.1) WRITE(7,6350) WO,HO,XO,DO,RHO(2,3) GOTO 39 C *** COMPREHENSIVE INPUT FOR GENERAL ROOM; C *** NON-RECTANGULAR ROOMS, NON-RECTANGULAR SURFACES, OBSTRUCTIONS, ETC. C C DATA FOR ROOM C 31 WRITE(7,6005) C NUMBER OF DIFFERENT TYPES OF SURFACES READ(5,*) NC,NSW,ND,NWC,NWO,NWS,NSS WRITE(7,6000)NC,NSW,ND,NWC,NWO,NWS,NSS N0=NC+NSW N01=NC+1 N1=N0+ND N11=N1+1 N2=N1+NWC NT=N2+NWO NT1=NT+1 NTS=NT+NWS NSP1=NSS+1 IF (NWS.LE.3) GOTO 5170 WRITE(7,8007) STOP 5170 IF (NSS.LE.2) GOTO 5180 WRITE(7,8008) STOP 5180 IF (N1*2+NWC+NWO+NWS.LE.30) GOTO 5200 WRITE(7,8010) STOP C 5200 DO 10 I=1,NT C INTERNAL SURFACE DATA READ(5,*) (X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ=1, 1 2),RHO(1,I),NM(I),NN(I) WRITE(7,6010)I,(X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ=1, 1 2),RHO(1,I) WRITE(7,6020)NM(I),NN(I) IF (NM(I)*NN(I).LE.400) GOTO 5205 WRITE(7,8011) I STOP 5205 IF (RHO(1,I).GE.0..AND.RHO(1,I).LE.1.) GOTO 5210 WRITE(7,8012) I STOP 5210 DO 11 IZ=1,2 BETA(1,I,IZ)=RAD*BETA(1,I,IZ) 11 PSI(1,I,IZ)=RAD*PSI(1,I,IZ) C C NUMBER AND NAMES OF CUT-OUTS IF(I.LE.N2) GO TO 13 READ(5,*) ICO(I),(JCO(I,M),M=1,ICO(I)) IF (ICO(I).LT.1) GOTO 13 DO 5020 M=1,ICO(I) IF (JCO(I,M).GE.1.AND.JCO(I,M).LE.NTS) GOTO 5020 WRITE(7,8014) I STOP 5020 CONTINUE C C NUMBER AND NAMES OF OBSTRUCTORS C 13 READ(5,*) (IOB(1,I,J),(JOB(I,J,K),K=1,IOB(1,I,J)),J=1,NT) IDS(I)=1 IF(I.GT.N1) GO TO 10 C WINDOW DATA READ(5,*) IS,AK2(I),TAU(I),DW(I),IOH(I) IF (IS.LE.NSS) GOTO 5220 WRITE(7,8020) I STOP 5220 IF (TAU(I).GE.0..AND.TAU(I).LE.1.) GOTO 5230 WRITE(7,8022) I STOP 5230 IF (AK2(I).GE.0..AND.AK2(I).LE.1.) GOTO 5240 WRITE(7,8024) I STOP 5240 IF(IOH(I).EQ.0) GOTO 16 C OVERHANG DATA READ(5,*) OHW(I),RING(I),OHREF(I),KOH(I),OHTAU(I) IF (KOH(I).EQ.0) OHTAU(I)=0. IF (OHREF(I).GE.0..AND.OHREF(I).LE.1.) GOTO 5250 WRITE(7,8034) I STOP 5250 IF (OHTAU(I).GE.0..AND.OHTAU(I).LE.1.) GOTO 16 WRITE(7,8036) I STOP 16 IDS(I)=IS+1-1/(1+IS) ALPHA(I)=1. IF(I.GT.N0)ALPHA(I)=0. ISD(I)=1 WRITE(7,6030)IS,AK2(I),TAU(I),DW(I) IF(I.LE.NC.OR.I.GT.N0) GO TO 18 C SHEER CURTAIN READ(5,*) ALPHA(I) WRITE(7,6064) ALPHA(I) 18 IF(IOH(I).EQ.0) GOTO 10 IF(IOH(I).EQ.1) WRITE(7,6031) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.2) WRITE(7,6032) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.3) WRITE(7,6033) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.4) WRITE(7,6034) OHW(I),RING(I),OHREF(I) IF(IOH(I).EQ.5) WRITE(7,6035) OHW(I),RING(I),OHREF(I) IF(KOH(I).EQ.0) WRITE(7,6036) IF(KOH(I).EQ.1) WRITE(7,6037) OHTAU(I) IF(KOH(I).EQ.2) WRITE(7,6038) OHTAU(I) 10 CONTINUE C C DATA FOR WORKING SURFACES C DO 35 I=NT1,NTS READ(5,*) (X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ 1 =1,2),NM(I),NN(I) WRITE(7,6070)I,(X0(1,I,IZ),IZ=1,7),(BETA(1,I,IZ),PSI(1,I,IZ),IZ=1 1 ,2),NM(I),NN(I) IF (NM(I)*NN(I).LE.400) GOTO 5270 WRITE(7,8011) I STOP 5270 RHO(1,I)=1. DO 36 IZ=1,2 BETA(1,I,IZ)=RAD*BETA(1,I,IZ) 36 PSI(1,I,IZ)=RAD*PSI(1,I,IZ) READ(5,*) (IOB(1,I,J),(JOB(I,J,K),K=1,IOB(1,I,J)),J=1,NT) 35 CONTINUE C C DATA FOR OUTSIDE ENCLOSURES C NW(1)=0 IF(NSS.EQ.0)GO TO 39 C SCAN OVER ALL OUTSIDE ENCLOSURES DO 20 K=2,NSP1 C NUMBER OF SOLID WALLS IN OUTSIDE ENCLOSURE READ(5,*) NW(K) K2=K-1 WRITE(7,6040)K2,NW(K) NTK=NW(K) DO 20 I=1,NTK C SURFACE DATA READ (5,*) (X0(K,I,IZ),IZ=1,7),(BETA(K,I,IZ),PSI(K,I,IZ),IZ=1, 1 2),RHO(K,I) WRITE(7,6010)I,(X0(K,I,IZ),IZ=1,7),(BETA(K,I,IZ),PSI(K,I,IZ),IZ=1, 1 2),RHO(K,I) IF (RHO(K,I).GE.0..AND.RHO(K,I).LE.1.) GOTO 5290 WRITE(7,8005) I,K STOP 5290 DO 17 IZ=1,2 BETA(K,I,IZ)=RAD*BETA(K,I,IZ) 17 PSI(K,I,IZ)=RAD*PSI(K,I,IZ) C NUMBER OF OBSTRUCTORS 20 READ(5,*) (IOB(K,I,J),(KJOB(K,I,J),L=1,IOB(K,I,J)),J=1,NTK) C NUMBER OF OBSTRUCTORS BETWEEN INSIDE AND OUTSIDE SURFACES C DO 30 K=2,NSP1 NTK=NW(K) DO 30 I=1,NTS 30 READ(5,*) (KOB(K,I,J),(KKOB(K,I,L),L=1,KOB(K,I,J)),J=1,NTK) 39 WSOR=PSI(1,NT1,1) NW(1)=0 C C ILLUMINATION DATA C MFIRST=1 MLAST=1 MSTEP=1 NTSTEP=1 IDAY=1 GOTO(5400,34,5405)IWMOD C SPECIFIC SOLAR DATA ARE SUPPLIED 5400 READ(5,*) SOLC1,ETSOL,SKYC1,ETSKY,ICLD,RHOGR WRITE(7,6050)SOLC1,ETSOL,SKYC1,ETSKY,ICLD,RHOGR SKYC=ETSKY*SKYC1 IF(ICLD.EQ.1.OR.ICLD.EQ.3) GO TO 32 C CLEAR DAY READ(5,*) HSUND,PSID GOTO 5410 C GEOGRAPHICAL DATA ARE SUPPLIED; SCANS THROUGH MONTHS, HOURS 34 READ(5,*) ALAT,ALONG,TZN,ALT,ICLD,RHOGR READ(5,*) MFIRST,MLAST,MSTEP,IDAY,TFIRST,TLAST,DELTAT WRITE(7,6100) ALAT,ALONG,TZN,ALT,ICLD,RHOGR WRITE(7,6110) MFIRST,MLAST,MSTEP,IDAY,TFIRST,TLAST,DELTAT IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 42 READ(5,*) (W(M),BW(M),M=MFIRST,MLAST,MSTEP) WRITE(7,6120) DO 6125 M=MFIRST,MLAST,MSTEP IF (W(M).GE.0..AND.W(M).LE.10.) GOTO 6124 WRITE(7,8030) M STOP 6124 IF (BW(M).GE.0..AND.BW(M).LE..5) GOTO 6125 WRITE(7,8032) M STOP 6125 WRITE(7,6121) M,W(M),BW(M) 42 NTSTEP=(TLAST-TFIRST)/DELTAT+1.1 GOTO 32 C C SUN POSITION( ALTITUDE AND AZIMUTH) DATA ARE SUPPLIED C 5405 READ(5,*) HSUND,PSID,ALT,ICLD,RHOGR WRITE(7,6105) HSUND,PSID,ALT,ICLD,RHOGR IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 5410 READ(5,*) W(1),BW(1) 5410 HSUN=HSUND*RAD PSIS=PSID*RAD IISUN=1 32 IF(INID.EQ.1.AND.IOPP.NE.-1)RHO(2,2)=RHOGR IF (RHOGR.GE.0..AND.RHOGR.LE.1.) GOTO 5300 WRITE (6,8026) STOP C C EVALUATION OF DIRECTION COSINES C 5300 DO 40 K=1,NSP1 NTK=NTS IF(K.GE.2) NTK=NW(K) DO 40 I=1,NTK SBX=SIN(BETA(K,I,1)) CBX=COS(BETA(K,I,1)) SBY=SIN(BETA(K,I,2)) CBY=COS(BETA(K,I,2)) SPX=SIN(PSI(K,I,1)) CPX=COS(PSI(K,I,1)) SPY=SIN(PSI(K,I,2)) CPY=COS(PSI(K,I,2)) LL(K,I,1)=SBX*CPX LL(K,I,2)=SBX*SPX LL(K,I,3)=CBX LL(K,I,4)=SBY*CPY LL(K,I,5)=SBY*SPY LL(K,I,6)=CBY LL(K,I,7)=SBX*SPX*CBY-CBX*SBY*SPY LL(K,I,8)=CBX*SBY*CPY-SBX*CPX*CBY LL(K,I,9)=SBX*SBY*(SPY*CPX-CPY*SPX) 40 CONTINUE C C CALCULATION OF NODAL COORDINATES C DO 1100 I=1,NTS K=0 NM1=NM(I) NN1=NN(I) AA=X0(1,I,7)/(NM1*NN1) BB=X0(1,I,4)-X0(1,I,6)+X0(1,I,5) DO 1000 N=1,NN1 ETA=(NN1+.5-N)/NN1 CC=X0(1,I,4)-BB*ETA DO 1000 M=1,NM1 XSI=(M-.5)/NM1 XIP=X0(1,I,5)*ETA+CC*XSI YIP=X0(1,I,7)*ETA DO 1010 L=1,3 1010 XN(I,K+1,L)=X0(1,I,L)+LL(1,I,L)*XIP+LL(1,I,L+3)*YIP C CHECK WHETHER A NODAL POINT FALLS ONTO A CUT-OUT IF(I.LE.N2.OR.I.GT.NT) GOTO 1020 IICO=ICO(I) DO 1030 IC=1,IICO J=JCO(I,IC) XJP=0. YJP=0. DO 1040 L=1,3 XYP=XN(I,K+1,L)-X0(1,J,L) XJP=XJP+XYP*LL(1,J,L) 1040 YJP=YJP+XYP*LL(1,J,L+3) ETAJ=YJP/X0(1,J,7) IF(ABS(ETAJ-.5).GT..5) GOTO 1030 XSIJ=(XJP-X0(1,J,5)*ETAJ)/(X0(1,J,4)-(X0(1,J,4)-X0(1,J,6)+ 1 X0(1,J,5))*ETAJ) IF(ABS(XSIJ-.5).GT..5) GOTO 1030 GOTO 1000 1030 CONTINUE 1020 K=K+1 AIK(I,K)=AA*CC 1000 CONTINUE C NUMBER OF NODES ON SURFACE 1100 NUM(I)=K C C CALCULATION OF COORDINATES FOR OUTSIDE SURFACE OF WINDOW OPENING C (STORE THESE AUXILIARY SURFACES BEYOND WORKING SURFACE, I.E.@ NTS+I) C DO 45 I=1,N1 DO 43 J=1,3 X0(1,I+NTS,J)=X0(1,I,J)-(DW(I)+1.E-3)*LL(1,I,J+6) DO 43 K=1,3 JJ=J+3*(K-1) 43 LL(1,I+NTS,JJ)=LL(1,I,JJ) DO 45 J=4,7 45 X0(1,I+NTS,J)=X0(1,I,J) C C CALCULATION OF OUTSIDE SURFACE AREAS C IF(NSS.EQ.0) GOTO 49 DO 44 K=2,NSP1 NTK=NW(K) DO 44 I=1,NTK LUM(K,I)=0 44 A(K,I)=(X0(K,I,4)+X0(K,I,6)-X0(K,I,5))*X0(K,I,7)/2. C C GEOMETRIC DESCRIPTION OF WINDOW OVERHANGS C (STORED BEYOND OUTSIDE ENCLOSURES, I.E. 1+NSS+1) C 49 NOH=NSP1+1 IOHK=0 DO 490 I=1,N1 IF(IOH(I).EQ.0) GOTO 490 IN=I+NTS IOHK=IOHK+1 ISTART(I)=IOHK X2Y=X0(1,IN,5)/X0(1,IN,7) RTX2Y=SQRT(1.+X2Y*X2Y) X31Y=(X0(1,IN,4)-X0(1,IN,6))/X0(1,IN,7) RTX31Y=SQRT(1.+X31Y*X31Y) XLT=X0(1,IN,5)-RING(I)*(RTX2Y-X2Y) YLT=X0(1,IN,7)+RING(I) XTR=X0(1,IN,6)+RING(I)*(RTX31Y-X31Y) YTR=YLT IF(IOH(I).EQ.2.OR.IOH(I).EQ.3) GOTO 425 C TOP OVERHANG DO 410 L=1,3 410 X0(NOH,IOHK,L)=X0(1,IN,L)+LL(1,I,L)*XLT+LL(1,I,L+3)*YLT X0(NOH,IOHK,4)=OHW(I) X0(NOH,IOHK,5)=0. X0(NOH,IOHK,6)=OHW(I) X0(NOH,IOHK,7)=XTR-XLT C DO 420 L=1,3 LL(NOH,IOHK,L)=-LL(1,I,L+6) LL(NOH,IOHK,L+3)=LL(1,I,L) 420 LL(NOH,IOHK,L+6)=-LL(1,I,L+3) RHO(NOH,IOHK)=OHREF(I) IF(IOH(I).LE.1) GOTO 490 C IOHK=IOHK+1 425 XBL=-RING(I)*(RTX2Y+X2Y) YBL=-RING(I) XRB=X0(1,IN,4)+RING(I)*(RTX31Y+X31Y) YRB=YBL IF(IOH(I).EQ.3) GOTO 455 C LEFT AND RIGHT OVERHANGS DO 430 L=1,3 X0(NOH,IOHK,L)=X0(1,IN,L)+LL(1,I,L)*XBL+LL(1,I,L+3)*YBL 430 X0(NOH,IOHK+1,L)=X0(1,IN,L)+LL(1,I,L)*XTR+LL(1,I,L+3)*YTR X0(NOH,IOHK,4)=OHW(I) X0(NOH,IOHK,5)=0. X0(NOH,IOHK,6)=OHW(I) X0(NOH,IOHK,7)=(X0(1,IN,7)+2.*RING(I))*RTX2Y X0(NOH,IOHK+1,4)=OHW(I) X0(NOH,IOHK+1,5)=0. X0(NOH,IOHK+1,6)=OHW(I) X0(NOH,IOHK+1,7)=(X0(1,IN,7)+2.*RING(I))*RTX31Y C DO 440 L=1,3 LL(NOH,IOHK,L)=-LL(1,I,L+6) LL(NOH,IOHK,L+3)=(LL(1,I,L)*X2Y+LL(1,I,L+3))/RTX2Y LL(NOH,IOHK+1,L)=-LL(1,I,L+6) 440 LL(NOH,IOHK+1,L+3)=(LL(1,I,L)*X31Y-LL(1,I,L+3))/RTX31Y DO 450 L=1,3 L1=L+1-3*(L/3) L2=L+2-3*(L/2) LL(NOH,IOHK,L+6)=LL(NOH,IOHK,L1)*LL(NOH,IOHK,L2+3) 1 -LL(NOH,IOHK,L2)*LL(NOH,IOHK,L1+3) 450 LL(NOH,IOHK+1,L+6)=LL(NOH,IOHK+1,L1)*LL(NOH,IOHK+1,L2+3) 1 -LL(NOH,IOHK+1,L2)*LL(NOH,IOHK+1,L1+3) RHO(NOH,IOHK)=OHREF(I) IOHK=IOHK+1 RHO(NOH,IOHK)=OHREF(I) IF(IOH(I).EQ.2.OR.IOH(I).EQ.4) GOTO 490 C C BOTTOM OVERHANG IOHK=IOHK+1 455 DO 460 L=1,3 460 X0(NOH,IOHK,L)=X0(1,IN,L)+LL(1,I,L)*XRB+LL(1,I,L+3)*YRB X0(NOH,IOHK,4)=OHW(I) X0(NOH,IOHK,5)=0. X0(NOH,IOHK,6)=OHW(I) X0(NOH,IOHK,7)=XRB-XBL C DO 470 L=1,3 LL(NOH,IOHK,L)=-LL(1,I,L+6) LL(NOH,IOHK,L+3)=-LL(1,I,L) 470 LL(NOH,IOHK,L+6)=LL(1,I,L+3) RHO(NOH,IOHK)=OHREF(I) 490 CONTINUE ISTART(N1+1)=IOHK+1 ZENL0=0 C C START CYCLE FOR MONTHS AND HOURS (IF IWMOD=2) C DO 3400 MONTH=MFIRST,MLAST,MSTEP DO 3400 ITIME=1,NTSTEP IF(IWMOD.NE.2) GOTO 494 TIME=TFIRST + DELTAT*(ITIME-1) WRITE(7,6130) MONTH,TIME C C ZENITH LUMINANCE C 494 DO 495 K=2,NSP1 NK=NW(K) DO 495 I=1,NK 495 DSOL(K,I)=0. DO 496 I=1,NWS KMAX=NUM(I+NT) DO 496 K=1,KMAX 496 WSOL(I,K)=0. DO 497 I=1,NTS KMAX=NUM(I) DO 497 K=1,KMAX 497 LUMI(I,K)=0. GO TO(50,60,50,60),ICLD C ICLD=1.. OVERCAST SKY, SKY LUMINANCE DISTRIBUTION FROM C.I.E. C ICLD=2.. CLEAR SKY WITH DIRECT SUN INCLUDED. C ICLD=3.. UNIFORM-LUMIN@ANCE SKY, ZENITH VALUE FROM OVERCAST C.I.E. C ICLD=4.. CLEAR SKY ONLY 50 AK=5.*(3.-2.*RHOGR)/7.*(1/ICLD) AUX(1)=1./(1.+AK) AUX(2)=AUX(1)*AK SOLC=0. GOTO(57,58,59)IWMOD C LUMINANCE FROM INSOLATION DATA 57 ZENL=(1.+AK)*SKYC/((1.+.667*AK)) GRNDL=SKYC*RHOGR/ZENL GO TO 100 C LUMINANCE FROM GEOGRAPHICAL DATA C DETERMINE HSUN=SUN ALTITUDE ANGLE 58 CALL SOLAR(ALAT,ALONG,TZN,MONTH,TIME,IISUN,HSUN,PSIS) 59 ZENL=92.94*(.123+8.6*SIN(HSUN))*PI SKYC=(1.+.667*AK)/(1.+AK)*ZENL GRNDL=SKYC*RHOGR/ZENL IF(MONTH-MFIRST+ITIME.EQ.1) GOTO 100 C C SET UP LUMINANCES FOR OTHER TIMES IF NOT CLEAR SKY C (FOR NON-CLEAR SKY CALCUL@ATIONS NEED TO BE DONE ONLY ONCE, AS C GEOMETRY AND SKY DISTRIBUTION REMAIN UNCHANGED; ONLY THE VALUE C OF ZENL=ZENITH LUMINANCE CHANGES) C C RATIO OF NEW AND OLD ZENITH LUMINANCES ZR=ZENL/ZENL0 ZENL0=ZENL C OUTSIDE LUMINANCES IF (NSS.LE.0) GOTO 2405 DO 2400 K=1,NSS NKT=NW(K+1) DO 2400 I=1,NKT 2400 LUM(K+1,I)=ZR*LUM(K+1,I) 2405 IF(IOHK.EQ.0) GOTO 2415 C OVERHANG LUMINANCES DO 2410 II=1,IOHK 2410 LUM(NOH,II)=ZR*LUM(NOH,II) C INSIDE AND WORKING SURFACES 2415 DO 2420 I=1,NTS KMAX=NUM(I) DO 2420 K=1,KMAX 2420 LUMF(I,K)=ZR*LUMF(I,K) GOTO 1405 C C ESTABLISH SUN POSITION FOR THE GEOGRAPHICAL DATA GIVEN 60 IF(IWMOD.EQ.2)CALL SOLAR(ALAT,ALONG,TZN,MONTH,TIME, 1 IISUN,HSUN,PSIS) BETAS=PI/2.-HSUN HSUND=HSUN/RAD PSID=PSIS/RAD IF(IWMOD.EQ.2)WRITE(7,6140) HSUND,PSID C INTEGRATE NORMALIZED LUMINANCE OVER SKY P=SKYSUM(BETAS) AUX(3)=.27385*(.91+10*EXP(-3*BETAS)+.45*COS(BETAS)**2) IF (IWMOD.GT.1) GOTO 67 C IF SOLAR DATA KNOWN CALCULATE ZENL FROM SKY INTEGRAL SOLC=ETSOL*SOLC1 ES=SOLC/COS(BETAS) ZENL=SKYC/P TERM=SOLC/(SOLC+SKYC) C IF TOO LITTLE DIFFUSE LIGHT, FAULTY INPUT DATA HAVE BEEN GIVEN IF(TERM.LT..9) GOTO 33 WRITE(7,8000) STOP C ROUTINE TO CALCULATE ZENL AS FUNCTION OF TURBIDITY, VAPOR & SOL.ALT. 67 CALL ZENLCD(P,W(MONTH),BW(MONTH),HSUN,ZENL) C ROUTINE TO CALCULATE DIRECT SUN FROM SIMILAR DATA IF(ICLD.EQ.2)SOLC=DIR(MONTH,IDAY,BW(MONTH),W(MONTH),ALT,HSUN) IF(ICLD.EQ.4)SOLC=0. C USE SKY INTEGRAL AND ZENL TO GET DIFFUSE SKY ILLUMINATION SKYC=P*ZENL 33 II=0 ES=SOLC/COS(BETAS) GRNDL=(SKYC+SOLC)*RHOGR/ZENL IF(NSS.LE.0) GOTO 70 C INITIALIZATION OF SUNNY-NODES IDENTIFIERS FOR OUTSIDE SURFACES DO 38 K=2,NSP1 NK=NW(K) DO 38 J=1,NK II=II+1 KSUN(K,J)=0 DSOL(K,J)=0. DO 37 M=1,11 DO 37 N=1,11 KD(II,M,N)=0. 37 CONTINUE 38 CONTINUE C C DIRECTION COSINES FOR VECTOR TOWARDS SUN C 70 LS(1)=SIN(BETAS)*COS(PSIS) LS(2)=SIN(BETAS)*SIN(PSIS) LS(3)=COS(BETAS) C C CALCULATION OF DIRECT SUNSHINE ILLUMINATION INSIDE ROOM C WRITE(7,6067) DO 1270 IW=1,N1 ISUN(IW)=0 CBW=LS(1)*LL(1,IW,7)+LS(2)*LL(1,IW,8)+LS(3)*LL(1,IW,9) C DOES SUN HIT WINDOW? IF(CBW.GE.-1.E-4) GOTO 1270 C ANY OBSTRUCTION BETWEEN WINDOW AND SUN? KK=IDS(IW) IF(KK.LE.0) GOTO 1240 NWK=NW(KK) DO 1250 L=1,3 1250 X(L)=X0(1,IW,L) DO 1230 J=1,NWK IF(IOB(KK,1,J).LT.0) GOTO 1230 CALL SECT(X,LS,KK,J,KUT) IF(KUT.EQ.1) GOTO 1270 1230 CONTINUE 1240 ISUN(IW)=1 WRITE(7,6068) IW,ALPHA(IW) 1270 CONTINUE C IF(N0.EQ.0) GOTO 1300 DO 1200 I=1,NTS C COULD SUN POSSIBLY SHINE ON THIS SURFACE? CBS=LS(1)*LL(1,I,7)+LS(2)*LL(1,I,8)+LS(3)*LL(1,I,9) IF(CBS.LE.1.E-4) GOTO 1200 KMAX=NUM(I) DO 1210 IW=1,N0 IF(ISUN(IW).EQ.0) GOTO 1210 IF(IOB(1,I,IW).LT.0) GOTO 1210 CBW=LS(1)*LL(1,IW,7)+LS(2)*LL(1,IW,8)+LS(3)*LL(1,IW,9) C ATTENUATED AMOUNT OF DIRECT SUN LIGHT TO BE ADDED TO LUMINANCE WS=RHO(1,I)*ES*CBS*ALPHA(IW)*TAUB(IW,-CBW) DO 1260 K=1,KMAX DO 1220 L=1,3 1220 X(L)=XN(I,K,L) CALL SECT(X,LS,1,IW,KUT) C DOES RAY FROM NODE TO SUN PASS THROUGH WINDOW? IF(KUT.NE.1) GOTO 1260 C BACK SIDE OF WINDOW ALSO? CALL SECT(X,LS,1,IW+NTS,KUT) IF(KUT.NE.1) GOTO 1260 C HOW ABOUT INTERNAL OBSTRUCTIONS? IF(IOB(1,I,IW).EQ.0) GOTO 1255 IIOB=IOB(1,I,IW) DO 1245 IO=1,IIOB JO=JOB(I,IW,IO) CALL SECT(X,LS,1,JO,KUT) IF(KUT.EQ.1) GOTO 1260 1245 CONTINUE C ANY OVERHANG IN THE WAY? 1255 IF(IOH(IW).EQ.0) GOTO 1257 II1=ISTART(IW) II2=ISTART(IW+1)-1 DO 1256 II=II1,II2 CALL SECT(X,LS,NOH,II,KUT) IF(KUT.NE.1) GOTO 1256 C IF SO, IS OVERHANG SEMI-TRANSPARENT? (IN THAT CASE ATTENUATE MORE) IF(KOH(IW).NE.1) GOTO 1260 WS=WS*OHTAU(IW) GOTO 1257 1256 CONTINUE C ANY EXTERNAL OBSTRUCTION IN THE WAY? 1257 KK=IDS(IW) IF (KK.LE.0) GOTO 1275 NWK=NW(KK) DO 1258 J=1,NWK IF (KOB(KK,I,J).LT.0) GOTO 1258 CALL SECT(X,LS,KK,J,KUT) IF (KUT.EQ.1) GOTO 1260 1258 CONTINUE C DIRECT SUN IS INITIAL VALUE OF LUMINANCE 1275 LUMI(I,K)=WS 1260 CONTINUE 1210 CONTINUE 1200 CONTINUE C STORE THESE VALUES ALSO IN WSOL FOR PLOTTING ROUTINE DO 1280 I=NT1,NTS KMAX=NUM(I) DO 1280 K=1,KMAX 1280 WSOL(I-NT,K)=LUMI(I,K) C C CALCULATION OF DIRECT SUNSHINE ON CURTAINED AND C DIFFUSE WINDOWS C IF(N1.EQ.NC) GOTO 1350 1300 DO 1310 I=N01,N1 C DOES WINDOW GET SUNSHINE? IF(ISUN(I).EQ.0) GOTO 1310 CBW=LS(1)*LL(1,I,7)+LS(2)*LL(1,I,8)+LS(3)*LL(1,I,9) C CALCULATE ATTENUATED CONTRIBUTION TO LUMINANCE (NOTE CBW<0) WS=-ES*CBW*(1.-ALPHA(I))*TAUB(I,-CBW) KMAX=NUM(I) C ANY OVERHANG IN THE WAY FOR THIS NODE? DO 1320 K=1,KMAX TT=1. IF(IOH(I).EQ.0) GOTO 1320 DO 1330 L=1,3 1330 X(L)=XN(I,K,L) II1=ISTART(I) II2=ISTART(I+1)-1 DO 1340 II=II1,II2 CALL SECT(X,LS,NOH,II,KUT) IF(KUT.NE.1) GOTO 1340 TT=0. C IS OVERHANG SEMI-TRANSPARENT? IF(KOH(I).EQ.1) TT=OHTAU(I) GOTO 1320 1340 CONTINUE 1320 LUMI(I,K)=WS*TT 1310 CONTINUE C C CALCULATION OF DIRECT SUNSHINE ILLUMINATION FOR OUTSIDE ENCLOSURE SURFACES C 1350 IF(NSS.LE.0)GO TO 210 II=0 DO 90 K=2,NSP1 NK=NW(K) DO 90 J=1,NK II=II+1 FR(J)=0. C ARE SUNNY AREAS POSSIBLE ON THIS SURFACE? CBJ=LS(1)*LL(K,J,7)+LS(2)*LL(K,J,8)+LS(3)*LL(K,J,9) IF(CBJ.LE.1.E-4) GO TO 90 C C BREAK UP SURFACE INTO 11 BY 11 NODES AND CHECK FOR SUN (SEND RAYS TO SUN) DO 95 M=1,11 ETAM=(M-1.)/10. DO 95 N=1,11 ETAN=(N-1.)/10. XPP=X0(K,J,4)-(X0(K,J,4)-X0(K,J,6)+X0(K,J,5))*ETAM XP=X0(K,J,5)*ETAM+XPP*ETAN YP=X0(K,J,7)*ETAM DO 94 L=1,3 94 X(L)=X0(K,J,L)+LL(K,J,L)*XP+LL(K,J,L+3)*YP C ANYTHING BLOCKING THE SUN? DO 93 JJ=1,NK IF(J.EQ.JJ) GO TO 93 CBL=LS(1)*LL(K,JJ,7)+LS(2)*LL(K,JJ,8)+LS(3)*LL(K,JJ,9) IF(CBL.GE.-1.E-4)GO TO 93 CALL SECT (X,LS,K,JJ,KUT) IF(KUT.EQ.1) GO TO 95 93 CONTINUE KD(II,M,N)=1 C SUNNY FRACTION OF SURFACE FR(J)=FR(J)+C(M,11)*C(N,11)*XPP 95 CONTINUE C IF(FR(J).LE.1.E-4)GO TO 90 FR(J)=FR(J)*2/(X0(K,J,4)+X0(K,J,6)-X0(K,J,5)) FS(II)=(1.-FR(J))/FR(J) KSUN(K,J)=1 C CONTRIBUTION OF DIRECT SUN TO OUTSIDE SURFACE LUMINANCE DSOL(K,J)=FR(J)*CBJ*ES*RHO(K,J) 90 CONTINUE C C CALCULATION OF F I-J IN OUTSIDE ENCLOSURES C 100 ZENL0=ZENL IF(NSS.LE.0) GOTO 210 DO 400 K=2,NSP1 NMW=NW(K) CALL CONFAC(K,NMW,F) 400 CONTINUE C C LUMINANCES FOR EXTERNAL REFLECTIONS (SOLUTION OF SIMULTANEOUS EQS. C BY MATRIX INVERSION) C DO 220 K=2,NSP1 NKT=NW(K) DO 190 I=1,NKT DO 180 J=1,NKT DIJ=(I/J)*(J/I) 180 COEF(I,J)=DIJ/RHO(K,J)-F(K,I,J) 190 COFB(I,1)=ZENL*F(K,I,NKT+1)+DSOL(K,I)/RHO(K,I) IF (NKT.EQ.1) GOTO 230 CALL MATINV(COEF,NKT,COFB,1,DETERM) DO 200 I=1,NKT 200 LUM(K,I)=COFB(I,1) GOTO 220 230 LUM(K,1)=COFB(1,1)/COEF(1,1) 220 CONTINUE C C CALCULATION OF OVERHANG LUMINANCES C (USING MONTE CARLO, ASSUMING NO INTERFERENCE BETWEEN OVERHANGS) C 210 NBOH=NBUNDL/5 DO 1495 IW=1,N1 IF(IOH(IW).EQ.0) GOTO 1495 K=IDS(IW) C NUMBER OF OVERHANG SURFACES FOR THIS WINDOW II1=ISTART(IW) II2=ISTART(IW+1)-1 DO 1450 II=II1,II2 SUM=0. C SEND OUT BUNDLES DO 1455 NBUN=1,NBOH C ************************ RANDOM NUMBER GENERATOR ********************** DO 112,INDEX=1,4 112 RS(INDEX)=RANDOM(ISEED) C CHOOSE DIRECTION SINB=SQRT(RS(1)) COSB=SQRT(1.-RS(1)) THET=2.*PI*RS(2) COST=COS(THET) SINT=SIN(THET) C CHOOSE EMISSION POINT DO 1460 L=1,3 X(L)=X0(NOH,II,L)+LL(NOH,II,L)*X0(NOH,II,4)*RS(3) 1 +LL(NOH,II,L+3)*X0(NOH,II,7)*RS(4) Y(L)=(LL(NOH,II,L)*COST+LL(NOH,II,L+3)*SINT)*SINB 1 +LL(NOH,II,L+6)*COSB 1460 YM(L)=-Y(L) IF(K.EQ.0) GOTO 1466 C DOES IT HIT OUTSIDE SURFACE OF WINDOW? NWK=NW(K) CALL SECT(X,YM,1,IW+NTS,KUT) IF(KUT.NE.1) GOTO 1462 SUM=SUM+RHO(1,IW)/RHO(K,1)*LUM(K,1) GOTO 1455 C ANY OUTSIDE SURFACE? 1462 IF (NSS.LT.1) GOTO 1466 DO 1465 IR=1,NWK CALL SECT(X,Y,K,IR,KUT) IF(KUT.NE.1) GOTO 1465 SUM=SUM+LUM(K,IR) GOTO 1455 1465 CONTINUE C HITS THE SKY, (OR GROUND IF Y POINTS BELOW HORIZ) 1466 SUM=SUM+RLPZ(ICLD,Y)*ZENL 1455 CONTINUE LUM(NOH,II)=OHREF(IW)*SUM/NBOH IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 1470 C DOES DIRECT SUN HIT INSIDE OVERHANG SURFACES? IF(ISUN(IW).EQ.0) GOTO 1470 CBS=LS(1)*LL(NOH,II,7)+LS(2)*LL(NOH,II,8)+LS(3)*LL(NOH,II,9) IF(CBS.LT.1.E-4) GOTO 1470 IF(K.EQ.0) GOTO 1476 DO 1475 IR=1,NWK CALL SECT(X,LS,K,IR,KUT) IF(KUT.EQ.1) GOTO 1470 1475 CONTINUE 1476 LUM(NOH,II)=LUM(NOH,II)+OHREF(IW)*ES*CBS 1470 IF(KOH(IW).NE.2) GOTO 1450 C FOR TRANSLUCENT OVERHANG ALSO NEED LUMINANCE CONTRIBUTION FROM OUTSIDE C OF OVERHANG SURFACE; ROUTINE SAME AS THE INSIDE (SEE ABOVE) SUM=0. DO 1480 NBUN=1,NBOH C ************************ RANDOM NUMBER GENERATOR ********************** DO 113 INDEX=1,4 113 RS(INDEX)=RANDOM(ISEED) SINB=SQRT(RS(1)) COSB=SQRT(1.-RS(1)) THET=2.*PI*RS(2) COST=COS(THET) SINT=SIN(THET) DO 1485 L=1,3 X(L)=X0(NOH,II,L)+LL(NOH,II,L)*X0(NOH,II,4)*RS(3) 1 +LL(NOH,II,L+3)*X0(NOH,II,7)*RS(4) 1485 Y(L)=-(LL(NOH,II,L)*COST+LL(NOH,II,L+3)*SINT)*SINB 1 -LL(NOH,II,L+6)*COSB IF(K.EQ.0) GOTO 1491 DO 1490 IR=1,NWK CALL SECT(X,Y,K,IR,KUT) IF(KUT.NE.1) GOTO 1490 SUM=SUM+LUM(K,IR) GOTO 1480 1490 CONTINUE 1491 SUM=SUM+RLPZ(ICLD,Y)*ZENL 1480 CONTINUE LUM(NOH,II)=LUM(NOH,II)+OHTAU(IW)*SUM/NBOH IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 1450 IF(ISUN(IW).EQ.0) GOTO 1450 CBS=-LS(1)*LL(NOH,II,7)-LS(2)*LL(NOH,II,8)-LS(3)*LL(NOH,II,9) IF(CBS.LT.1.E-4) GOTO 1450 IF(K.EQ.0) GOTO 1451 DO 1505 IR=1,NWK CALL SECT(X,LS,K,IR,KUT) IF(KUT.EQ.1) GOTO 1450 1505 CONTINUE 1451 LUM(NOH,II)=LUM(NOH,II)+OHTAU(IW)*ES*CBS 1450 CONTINUE 1495 CONTINUE C C OUTSIDE ILLUMINATION ON CURTAINED AND DIFFUSE WINDOWS C (USING MIXED METHOD.. INTEGRATION AND MONTE CARLO) C IF(NC.EQ.N1) GOTO 1600 DO 1500 I=N01,N1 IO=I+NTS XY1=(X0(1,IO,6)-X0(1,IO,5))/X0(1,IO,4) ID1=1 IF(ABS(XY1-1.).LT.1.E-3) ID1=0 AX=XY1*ID1 HI=0. K=IDS(I) C INTEGRATE OVER ALL DIRECTIONS DO 1550 IETA=1,21 ETA=(IETA-1.)/20. SETA=SQRT(1.-ETA*ETA) TT=C(IETA,21)*ETA*TAUB(I,ETA) DO 1550 IMU=1,21 AMU=(IMU-1.)/20. SWY=SETA*COS(2.*PI*AMU) SWZ=SETA*SIN(2.*PI*AMU) C EMISSION LOCATION FIXED BY RANDOM NUMBERS C ************************ RANDOM NUMBER GENERATOR ********************** DO 114,INDEX=1,4 114 RS(INDEX)=RANDOM(ISEED) YY=(1.-SQRT(1.-(1.-XY1*XY1)*RS(1)))/(1.-AX)+(1-ID1)*RS(1) XX=X0(1,IO,5)*YY+X0(1,IO,4)*(1.-(1.-XY1)*YY)*RS(2) DO 1520 L=1,3 X(L)=X0(1,IO,L)+LL(1,I,L)*XX+LL(1,I,L+3)*YY*X0(1,IO,7) 1520 Y(L)=-ETA*LL(1,I,L+6)+SWY*LL(1,I,L+3)+SWZ*LL(1,I,L) C ANY OVERHANG IN THE WAY? IF(IOH(I).EQ.0) GOTO 1545 II1=ISTART(I) II2=ISTART(I+1)-1 DO 1515 II=II1,II2 CALL SECT(X,Y,NOH,II,KUT) IF(KUT.NE.1) GOTO 1515 SS=LUM(NOH,II) GOTO 1550 1515 CONTINUE C ANY OUTSIDE SURFACE IN THE WAY? 1545 IF(K.EQ.0) GOTO 1530 NWK=NW(K) DO 1540 J=1,NWK IF(KOB(K,I,J)) 1540,1535,1525 1525 JR=KKOB(K,I,J) SS=LUM(K,JR) CALL SECT(X,Y,K,JR,KUT) IF(KUT.EQ.1) GOTO 1550 1535 CALL SECT(X,Y,K,J,KUT) IF(KUT.NE.1) GOTO 1540 SS=LUM(K,J) GOTO 1550 1540 CONTINUE C HITS THE SKY, (OR GROUND IF Y POINTS BELOW HORIZ) 1530 SS=ZENL*RLPZ(ICLD,Y) 1550 HI=HI+TT*SS*C(IMU,21) KMAX=NUM(I) DO 1500 K=1,KMAX 1500 LUMI(I,K)=LUMI(I,K)+2.*(1.-ALPHA(I)*(1/ISD(I)))*HI C C OUTSIDE ILLUMINATION ON INSIDE NODES C (INCLUDES INSIDE LUMINANCE OF WINDOWS, IN PARTICULAR DIFFUSING WINDOWS) C 1600 DO 1750 I=1,NTS KMAXI=NUM(I) FRAC=RHO(1,I)/PI DO 1750 J=1,N1 C CAN I SEE J? IF(IOB(1,I,J).LT.0) GOTO 1750 KMAXJ=NUM(J) K=IDS(J) IF(K.GT.0.AND.J.LE.NC) GOTO 1610 IF(K.GT.0.AND.ISD(J).EQ.1) GOTO 1610 C AUXILIARY DATA FOR COMPUTER EFFICIENCY, FOR CASE OF NO OUTSIDE ILLUMINATION C POSSIBLE ON THAT NODE K=4 LUM(4,1)=0. K1=30 K2=30 KOB(4,I,30)=-1 GOTO 1620 1610 K1=1 K2=NW(K) 1620 IF(IOB(1,I,J).GT.0) IIOB=IOB(1,I,J) C THE FOLLOWING ALGORITHMS ARE FOR CASES WITH & WITHOUT INTERNAL OBSTRUCTORS DO 1660 KI=1,KMAXI C START CALCULATIONS FOR NODE KI ON SURFACE I DO 1630 L=1,3 1630 X(L)=XN(I,KI,L) HIT=0. DO 1642 KJ=1,KMAXJ C DIRECTION, COSINES*S, AND DISTANCE-SQUARED TO NODE KJ ON J (WINDOW) SCB1=0. SCB2=0. SSQ=0. DO 1650 L=1,3 Y(L)=XN(J,KJ,L)-X(L) SCB1=SCB1+Y(L)*LL(1,I,L+6) SCB2=SCB2-Y(L)*LL(1,J,L+6) 1650 SSQ=SSQ+Y(L)*Y(L) IF(IOB(1,I,J).EQ.0) GOTO 1680 C ANY INTERNAL OBSTRUCTORS IN THE WAY? DO 1710 IO=1,IIOB JO=JOB(I,J,IO) CALL SECT(X,Y,1,JO,KUT) IF(KUT.EQ.1) GOTO 1642 1710 CONTINUE 1680 S=SQRT(SSQ) DO 1675 L=1,3 1675 Y(L)=Y(L)/S ICOS1=1.+SCB1/S ICOS2=1.+SCB2/S C CONFIGURATION FACTOR (CORRECTED BY FUDGE FACTOR FOR CLOSE NODES) FIJ=ICOS1*ICOS2*SCB1*SCB2/(SSQ*SSQ)*AIK(J,KJ) FIJ=FIJ/(1.+.6*FIJ*FIJ) C ASSUME OUTSIDE-WINDOW-WALL LUMINANCE FOR WINDOW FRAME HI=LUM(K,1) CB2=SCB2/S C DOES RAY GET TO THE OUTSIDE? CALL SECT(X,Y,1,J+NTS,KUT) SDM=1. IF(IOB(1,I,J).LE.0.AND.KUT.NE.1) GOTO 1640 IF(KUT.NE.1) GOTO 1635 C DOES RAY THROUGH WINDOW HIT ANY OUTSIDE SURFACE? DO 1670 JKL=K1,K2 JK=JKL IF(KOB(K,I,JK)) 1670,1625,1615 1615 JKK=KKOB(K,I,JK) CALL SECT(X,Y,K,JKK,KUT) IF(KUT.NE.1) GOTO 1625 JK=JKK GOTO 1605 1625 CALL SECT(X,Y,K,JK,KUT) IF(KUT.NE.1) GOTO 1670 1605 HI=LUM(K,JK) IF(ICLD.NE.2.OR.KSUN(K,JK).NE.1) GOTO 1635 C IS POINT ON OUTSIDE SURFACE IN SUN OR IN SHADE? C INTERPOLATE OUTSIDE LUMINANCES IF CLOSE TO SUN/SHADE DIVIDE II=JK KM1=K-1 DO 1655 K3=1,KM1 1655 II=II+NW(K3) TL=LUM(K,JK)-DSOL(K,JK) TH=LUM(K,JK)+FS(II)*DSOL(K,JK) ETAM=YDP/X0(K,JK,7) XPP=X0(K,JK,4)-(X0(K,JK,4)-X0(K,JK,6)+X0(K,JK,5))*ETAM ETAN=(XDP-X0(K,JK,5)*ETAM)/XPP TY=10*ETAM+1 LY=TY TX=10*ETAN+1 LX=TX F1=KD(II,LY,LX)+(TY-LY)*(KD(II,LY+1,LX)-KD(II,LY,LX)) F2=KD(II,LY,LX+1)+(TY-LY)*(KD(II,LY+1,LX+1)-KD(II,LY,LX+1)) PART=(F1+(TX-LX)*(F2-F1)) HI=TL+PART*(TH-TL) GOTO 1635 1670 CONTINUE C HITS SKY, (OR GROUND IF Y POINTS BELOW HORIZ) HI=ZENL*RLPZ(ICLD,Y) C ANY OVERHANG IN THE WAY? 1635 IF(IOH(J).EQ.0) GOTO 1640 II1=ISTART(J) II2=ISTART(J+1)-1 DO 1645 II=II1,II2 CALL SECT(X,Y,NOH,II,KUT) IF(KUT.NE.1) GOTO 1645 HI=LUM(NOH,II)+OHTAU(J)*HI*(2-KOH(J)) GOTO 1640 1645 CONTINUE C LUMINANCE CONTRIBUTION FROM OUTSIDE (HI) AND WINDOW ITSELF 1640 HIT=HIT+(HI*TAUB(J,CB2)*ALPHA(J)*(1/ISD(J))*SDM+LUMI(J,KJ))*FIJ 1642 CONTINUE 1660 LUMI(I,KI)=LUMI(I,KI)+FRAC*HIT 1750 CONTINUE C C LUMINANCES FOR INTERNAL REFLECTIONS C C KEEP DISTINCTION BETWEEN TOTAL (LUMF) AND EXTERNAL CONTRIBUTION (LUMI) DO 1760 I=1,NT KMAX=NUM(I) DO 1760 K=1,KMAX 1760 LUMF(I,K)=LUMI(I,K) C C GO THROUGH DESIRED NUMBER OF ITERATIONS DO 2200 ITER=1,NUMIT DO 1900 I=1,NT C IF WINDOW-REFLECTANCE IS SMALL, ITS CONTRIBUTION IS NEGLECTED C (FOR COMPUTER EFFICIENCY) C THE ALGORITHMS BELOW ARE BASICALLY THE SAME AS FOR OUTSIDE/WINDOW LUMINANCE C EFFECTS COMMENTED ON ABOVE IF(RHO(1,I).LE..15) GOTO 1900 FRAC=RHO(1,I)/PI KMAXI=NUM(I) DO 1805 KI=1,KMAXI 1805 DELF(KI)=0. DO 1890 J=N11,NT KMAXJ=NUM(J) IF(IOB(1,I,J)) 1890,1855,1850 C 1850 IIOB=IOB(1,I,J) 1855 DO 1880 KI=1,KMAXI DO 1860 L=1,3 1860 X(L)=XN(I,KI,L) DO 1880 KJ=1,KMAXJ SCB1=0. SCB2=0. SSQ=0. DO 1870 L=1,3 Y(L)=XN(J,KJ,L)-X(L) SCB1=SCB1+Y(L)*LL(1,I,L+6) SCB2=SCB2-Y(L)*LL(1,J,L+6) 1870 SSQ=SSQ+Y(L)*Y(L) IF(IOB(1,I,J).EQ.0) GOTO 1877 DO 1875 IO=1,IIOB CALL SECT(X,Y,1,JOB(I,J,IO),KUT) IF(KUT.EQ.1) GOTO 1880 1875 CONTINUE 1877 ICOS1=1.+SCB1/(1.+SSQ) ICOS2=1.+SCB2/(1.+SSQ) FIJ=SCB1*SCB2/(SSQ*SSQ)*AIK(J,KJ)*ICOS1*ICOS2 FIJ=FIJ/(1.+.6*FIJ*FIJ) DELF(KI)=DELF(KI)+FRAC*FIJ*LUMF(J,KJ) 1880 CONTINUE 1890 CONTINUE DO 1895 KI=1,KMAXI C ITERATION FINISHED FOR THIS SURFACE, IMPROVE VALUES FOR LUMF 1895 LUMF(I,KI)=LUMI(I,KI)+DELF(KI) 1900 CONTINUE 2200 CONTINUE C C CALCULATION OF ILLUMINATION AND DAYLIGHT FACTORS FOR WORKING SURFACE NODES C C ALGORITHMS ARE THE SAME AS FOR INTERNAL REFLECTIONS, BUT NO ITERATION, C BECAUSE OF THEIR IMPORTANCE, CONFIGURATION FACTORS FOR CLOSE NODES C ARE CALCULATED MORE ACCURATELY DO 1995 I=NT1,NTS KMAXI=NUM(I) DO 1905 KI=1,KMAXI 1905 DELF(KI)=0. DO 1990 J=N11,NT KMAXJ=NUM(J) IF(IOB(1,I,J)) 1990,1910,1950 1950 IIOB=IOB(1,I,J) 1910 DO 1980 KI=1,KMAXI DO 1960 L=1,3 1960 X(L)=XN(I,KI,L) DO 1980 KJ=1,KMAXJ SCB1=0. SCB2=0. SSQ=0. DO 1970 L=1,3 Y(L)=XN(J,KJ,L)-X(L) SCB1=SCB1+Y(L)*LL(1,I,L+6) SCB2=SCB2-Y(L)*LL(1,J,L+6) 1970 SSQ=SSQ+Y(L)*Y(L) IF(IOB(1,I,J).EQ.0) GOTO 1977 DO 1975 IO=1,IIOB CALL SECT(X,Y,1,JOB(I,J,IO),KUT) IF(KUT.EQ.1) GOTO 1980 1975 CONTINUE 1977 ICOS1=SCB1/(1.+SSQ)+1. ICOS2=SCB2/(1.+SSQ)+1. FIJ=ICOS1*ICOS2*SCB1*SCB2/(SSQ*SSQ)*AIK(J,KJ) IF(AIK(J,KJ)/SSQ.LT.1.) GOTO 2160 KKJ=KJ 2100 L1=(KKJ-1)/NM(J) L2=KKJ-L1*NM(J) IF(J.LE.N2) GOTO 2140 ETA=1.-(2.*L1+1.)/(2.*NN(J)) XSI=(2.*L2-1.)/(2.*NM(J)) CC=X0(1,J,4)-(X0(1,J,4)+X0(1,J,5)-X0(1,J,6))*ETA XJP=X0(1,J,5)*ETA+CC*XSI YJP=X0(1,J,7)*ETA XYZ=0. DO 2150 L=1,3 2150 XYZ=XYZ+X0(1,J,L)+LL(1,J,L)*XJP+LL(1,J,L+3)*YJP-XN(J,KJ,L) IF(ABS(XYZ).LT.1E-4) GOTO 2140 KKJ=KKJ+1 GOTO 2100 2140 FIJ=0. DO 2110 J1=1,3 ETA=1.-(2.*L1-(.5-J1)/1.5)/(2.*NN(J)) CC=X0(1,J,4)-(X0(1,J,4)+X0(1,J,5)-X0(1,J,6))*ETA AP=X0(1,J,7)*CC/(9.*NN(J)*NM(J)) DO 2110 J2=1,3 XSI=(2.*L2+(.5-J2)/1.5)/(2.*NM(J)) SCB1=0. SCB2=0. SSQ=0. XJP=X0(1,J,5)*ETA+CC*XSI YJP=X0(1,J,7)*ETA DO 2120 L=1,3 Y(L)=X0(1,J,L)+LL(1,J,L)*XJP+LL(1,J,L+3)*YJP-X(L) SCB1=SCB1+Y(L)*LL(1,I,L+6) SCB2=SCB2-Y(L)*LL(1,J,L+6) 2120 SSQ=SSQ+Y(L)*Y(L) ICOS1=SCB1/(1.+SSQ)+1. ICOS2=SCB2/(1.+SSQ)+1. 2110 FIJ=FIJ+ICOS1*ICOS2*SCB1*SCB2/(SSQ*SSQ)*AP 2160 DELF(KI)=DELF(KI)+FIJ*LUMF(J,KJ) 1980 CONTINUE 1990 CONTINUE C ADD THE INTERNAL-REFLECTION CONTRIBUTION TO THE WORKING SURFACE C SURFACE ILLUMINATION (WHICH SO FAR CONTAINS DIRECT SUN AND C CONTRIBUTION FROM OUTSIDE (INCLUDING SKY, EXTERNAL REFLECTIONS, C AND DIFFUSE WINDOW LIGHT, SUCH AS TRANSMISSION THROUGH CURTAINS, C SLATTED BLINDS, DIFFUSING WINDOWS) DO 1995 KI=1,KMAXI LUMF(I,KI)=LUMI(I,KI)+DELF(KI)/PI 1995 DAYF(I-NT,KI)=100.*(LUMF(I,KI)-WSOL(I-NT,KI))/SKYC C C PRINT WALL NODAL DATA AND LUMINANCES C 1405 WRITE(7,6250) SOLC,SKYC,ZENL IF(IWRITE.EQ.1) GOTO 1440 C C PRINTING OF LIGHT EXCHANGE FACTORS IN OUTSIDE ENCLOSURES C DO 1410 K=2,NSP1 NMW=NW(K) NM1=NMW+1 K4=K-1 WRITE(7,6560)K4 WRITE(7,6570)(J,J=1,NM1) WRITE(7,6572) DO 1410 I=1,NMW 1410 WRITE(7,6580) I,(F(K,I,J),J=1,NM1) C C PRINTING OF OUTSIDE LUMINANCES C DO 1420 K=1,NSS NKT=NW(K+1) WRITE(7,7010)K WRITE(7,7020)(I,I=1,NKT) 1420 WRITE(7,6200)(LUM(K+1,I),I=1,NKT) DO 1415 IW=1,N1 IF(IOH(IW).EQ.0) GOTO 1415 K=IDS(IW) II1=ISTART(IW) II2=ISTART(IW+1)-1 WRITE(7,6041) IW IF(IOH(IW).EQ.1) WRITE(7,6042)(LUM(NOH,II),II=II1,II2) IF(IOH(IW).EQ.2) WRITE(7,6043)(LUM(NOH,II),II=II1,II2) IF(IOH(IW).EQ.3) WRITE(7,6044)(LUM(NOH,II),II=II1,II2) IF(IOH(IW).EQ.4) WRITE(7,6045)(LUM(NOH,II),II=II1,II2) IF(IOH(IW).EQ.5) WRITE(7,6046)(LUM(NOH,II),II=II1,II2) 1415 CONTINUE WRITE(7,6710) DO 1400 I=1,NT KP=NM(I) KS=1+(KP-1)/12 L=NN(I) DO 1400 L1=1,L,KS KL=1+(L1-1)*KP KR=L1*KP IF(KL.GT.NUM(I)) GOTO 1440 IF(KR.GT.NUM(I)) KR=NUM(I) WRITE(7,6720)(K,K=KL,KR,KS) WRITE(7,6730)I,(XN(I,K,1),K=KL,KR,KS) WRITE(7,6740)(XN(I,K,2),K=KL,KR,KS) WRITE(7,6750)(XN(I,K,3),K=KL,KR,KS) WRITE(7,6760)(AIK(I,K),K=KL,KR,KS) 1400 WRITE(7,6770)(LUMF(I,K),K=KL,KR,KS) C C PRINT WORKING SURFACE NODAL DATA, ILLUMINATION, DAYLIGHT FACTOR C 1440 WRITE(7,6715) WRITE(7,6716) IF(INID.EQ.1)WRITE(7,6717) IF(INID.EQ.1)NT1=NWINDO+7 DO 1430 I=NT1,NTS KP=NM(I) KS=1+(KP-1)/12 L=NN(I) DO 1430 L1=1,L,KS KL=1+(L1-1)*KP KR=L1*KP WRITE(7,6720)(K,K=KL,KR,KS) IF(INID.EQ.1)WRITE(7,6730)I,(ABS(XN(I,K,1)),K=KL,KR,KS) IF(INID.EQ.2)WRITE(7,6730)I,(XN(I,K,1),K=KL,KR,KS) WRITE(7,6740)(XN(I,K,2),K=KL,KR,KS) WRITE(7,6750)(XN(I,K,3),K=KL,KR,KS) WRITE(7,6773) (LUMI(I,K),K=KL,KR,KS) WRITE(7,6774) (DELF(K)/PI,K=KL,KR,KS) WRITE(7,6775)(LUMF(I,K),K=KL,KR,KS) IF(INID.EQ.1)WRITE(7,6780)(DAYF(1,K),K=KL,KR,KS) 1430 IF(INID.EQ.2)WRITE(7,6780)(DAYF(I-NT,K),K=KL,KR,KS) IF(IPLOT.EQ.0) GOTO 3400 C C SET UP DATA FOR THE PLOTTING ROUTINE; CHOOSE LOCAL ORIGIN OF FIRST C WORKING SURFACE AS OVERALL ORIGIN FOR PLOT, AND STORE COORDINATES C FOR CORNERPOINTS@ (MUST BE RECTANGULAR) XY(1,1,1)=0. XY(1,1,2)=0. XY(1,2,1)=X0(1,NT1,5) XY(1,2,2)=X0(1,NT1,7) XY(1,3,1)=X0(1,NT1,6) XY(1,3,2)=X0(1,NT1,7) XY(1,4,1)=X0(1,NT1,4) XY(1,4,2)=0. C C SUBTRACT DIRECT SUN FROM ILLUMINATION FOR PLOTTING PURPOSES C AND RE-ORDER NODES KMAX=NUM(NT1) DO 3000 K=1,KMAX K1=(K-1)/NM(NT1) K2=K+(NN(NT1)-2*K1-1)*NM(NT1) 3000 LUMI(1,K)=LUMF(NT1,K2)-WSOL(1,K2) C C FIND ENDPOINTS FOR WINDOW PROJECTION ON FIRST WORKING SURFACE DO 2900 I=1,N1 DO 2800 L=1,3 PLTTMP(1,L)=X0(1,I,L) 2800 PLTTMP(2,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,6)+LL(1,I,L+3)*X0(1,I,7) DO 2900 I4=1,2 DO 2900 I2=1,2 XYW(I,I4,I2)=0. DO 2900 L=1,3 L2=L+3*(I2-1) 2900 XYW(I,I4,I2)=XYW(I,I4,I2)+(PLTTMP(I4,L)-X0(1,NT1,L))*LL(1,NT1,L2) C IF(NWS.EQ.1) GOTO 3200 C DO THE SAME AS ABOVE FOR ANY ADDITIONAL WORKING SURFACES C WITH RESPECT TO LOCAL AXES OF FIRST WORK SFC DO 3100 IW=2,NWS DO 3010 L=1,3 PLTTMP(1,L)=X0(1,NT+IW,L) PLTTMP(2,L)=PLTTMP(1,L)+LL(1,NT+IW,L)*X0(1,NT+IW,5)+ 1 LL(1,NT+IW,L+3)*X0(1,NT+IW,7) PLTTMP(3,L)=PLTTMP(1,L)+LL(1,NT+IW,L)*X0(1,NT+IW,6)+ 1 LL(1,NT+IW,L+3)*X0(1,NT+IW,7) 3010 PLTTMP(4,L)=PLTTMP(1,L)+LL(1,NT+IW,L)*X0(1,NT+IW,4) C DO 3020 I4=1,4 DO 3020 I2=1,2 XY(IW,I4,I2)=0. DO 3020 L1=1,3 LI2=L1+(I2-1)*3 3020 XY(IW,I4,I2)=XY(IW,I4,I2)+(PLTTMP(I4,L1)-X0(1,NT1,L1)) 1 *LL(1,NT1,LI2) C KMAX=NUM(NT+IW) DO 3100 K=1,KMAX K1=(K-1)/NM(NT+IW) K2=K+(NN(NT+IW)-2*K1-1)*NM(NT+IW) 3100 LUMI(IW,K)=LUMF(NT+IW,K2)-WSOL(IW,K2) C C THESE PLOTTING DATA ARE WRITTEN INTO A DATA FILE PLDAT WHICH C UPON JOB COMPLETION IS PERMANENTLY PLACED ON PSS-DISK C 3200 WRITE(8,8200) NT,NWS,ICLD,N0,N1,IPLOT WRITE(8,8200)(NM(I),NN(I),I=NT1,NTS) WRITE(8,8100) WSOR,SOLC,SKYC IF(ICLD.EQ.2.OR.ICLD.EQ.4) WRITE(8,8100) HSUND,PSID WRITE(8,8100) (((XYW(I,I4,I2),I2=1,2),I4=1,2),I=1,N1) DO 3150 I=1,NWS KMAX=NUM(NT+I) WRITE(8,8100)(LUMI(I,K),K=1,KMAX) 3150 WRITE(8,8100)((XY(I,I4,I2),I2=1,2),I4=1,4) IF(ICLD.EQ.1.OR.ICLD.EQ.3) GOTO 3400 C C ADDITIONAL PLOTTING DATA ARE GENERATED IF DIRECT SUN PENETRATES C INTO ENCLOSURE, VZ. THE IMAGE OF THE SUNNY AREA ON THE WORKING SURFACE C DO 3280 I=1,N0 DO 3205 IP=1,4 XSUN(I,IP,1)=0. 3205 XSUN(I,IP,2)=0. ANG(I)=0. IF(ISUN(I).EQ.0) GOTO 3280 DO 3210 L=1,3 PLTTMP(1,L)=X0(1,I,L) PLTTMP(2,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,5)+LL(1,I,L+3) 1 *X0(1,I,7) PLTTMP(3,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,6)+LL(1,I,L+3) 1 *X0(1,I,7) 3210 PLTTMP(4,L)=PLTTMP(1,L)+LL(1,I,L)*X0(1,I,4) IDENT=0 DO 3250 IP=1,4 DO 3220 L=1,3 X(L)=PLTTMP(IP,L) 3220 Y(L)=X(L)-(DW(I)+1E-3)*LL(1,I,L+6) IF(IDENT.NE.0) GOTO 3240 CALL SECT(X,LS,1,I+NTS,KUT) IF(KUT.EQ.1) IDENT=IP 3240 CALL SECT(X,LS,1,NT1,KUT) XSUN(I,IP,1)=XDP XSUN(I,IP,2)=YDP CALL SECT(Y,LS,1,NT1,KUT) 3250 CONTINUE IAP2=IDENT+2 IF (IAP2.GT.4) IAP2=IAP2-4 IIII=IDENT+1 IF (IIII.GT.4) IIII=IIII-4 JJJJ=IDENT+3 IF (JJJJ.GT.4) JJJJ=JJJJ-4 C DX=0. DY=0. DO 3260 L=1,3 DX=DX+LS(L)*LL(1,NT1,L) 3260 DY=DY+LS(L)*LL(1,NT1,L+3) ANG(I)=-180./PI*ATAN(DY/(DX+1E-8)) DO 3300 IP=1,ISNGRD DO 3300 L=1,ISNGRD 3300 IMG(I,IP,L)=0 ICNR(I,1)=IDENT ICNR(I,2)=IIII ICNR(I,3)=JJJJ ICNR(I,4)=IAP2 DO 3302 L=1,3 DO 3302 K=1,4 PLTTMP(K,L)=XSUN(I,ICNR(I,K),1)*LL(1,NT1,L) 2 +XSUN(I,ICNR(I,K),2)*LL(1,NT1,L+3) 2 +X0(1,NT1,L) 3302 CONTINUE DO 3390 M=1,ISNGRD DO 3390 N=1,ISNGRD DO 3304 L=1,3 X(L)=(PLTTMP(3,L)-PLTTMP(1,L))*(M-1)/(ISNGRD-1) 1 + PLTTMP(1,L) DUMMY=(PLTTMP(4,L)-PLTTMP(2,L))*(M-1)/(ISNGRD-1) 1 + PLTTMP(2,L) X(L)=(DUMMY-X(L))*(N-1)/(ISNGRD-1) + X(L) 3304 CONTINUE CALL SECT(X,LS,1,I+NTS,KUT) IF (KUT.NE.1) GOTO 3390 DO 3320 IWS=1,NWS IWSN=NT+IWS IIOB=IOB(1,IWSN,I) IF (IIOB.LT.1) GO TO 3320 DO 3310 K=1,IIOB JO=JOB(IWSN,I,K) CALL SECT(X,LS,1,JO,KUT) IF (KUT.NE.1) GO TO 3310 IF ((YDP+X0(1,JO,3)).LT.X0(1,IWSN,3)) GOTO 3310 IF (I.NE.4) GOTO 3390 GOTO 3390 3310 CONTINUE 3320 CONTINUE IF (IOH(I).EQ.0 .OR. KOH(I).EQ.1) GO TO 3370 IOHK1=ISTART(I) IOHK2=ISTART(I+1)-1 DO 3350 IOHK=IOHK1,IOHK2 CALL SECT(X,LS,NOH,IOHK,KUT) IF (KUT.EQ.1) GO TO 3390 3350 CONTINUE 3370 K=IDS(I) IF (K.LT.1) GO TO 3385 K2=NW(K) DO 3380 JK=1,K2 CALL SECT(X,LS,K,JK,KUT) IF (KUT.EQ.1) GO TO 3390 3380 CONTINUE 3385 IMG(I,M,N)=1 3390 CONTINUE 3280 CONTINUE WRITE(8,8100)(((XSUN(I,IP,J),J=1,2),IP=1,4),ANG(I),I=1,N0) WRITE(8,8200) (((IMG(I,M,N),N=1,ISNGRD),M=1,ISNGRD), 2 (ICNR(I,MM),MM=1,4),I=1,N0) C 3400 CONTINUE STOP C C6005 FORMAT(1H1,10X,19HINPUT DATA FOR ROOM/11X,19H*******************/) 6005 format(1H1,10X,'INPUT DATA FOR ROOM',/,11X,'*******************', 1/) 6000 FORMAT(/5X,33HNUMBER OF CLEAR WINDOWS ,I3/ 1 5X,33H SHEER CURTAIN WINDOWS ,I3/ 2 5X,33H DIFFUSE WINDOWS ,I3/ 3 5X,33H WALLS WITHOUT CUT-OUTS ,I3/ 4 5X,33H WALLS WITH CUT-OUTS ,I3/ 5 5X,33H WORKING SURFACES ,I3/ 6 5X,33H OUTSIDE ENCLOSURES ,I3///) 6010 FORMAT( 5X,7HSURFACE,I3/5X,10H----------// 1 10X,4HX0 =,F6.2,6H Y0 =,F6.2,6H Z0 =,F6.2/ 2 10X,4HX1 =,F6.2,6H X2 =,F6.2,6H X3 =,F6.2,5H Y =,F6.2/ 3 10X,7HBETAX =,F6.1,8H PSIX =,F6.1,9H BETAY =,F6.1, 4 8H PSIY =,F6.1/ 5 10X,14HREFLECTANCE =,F4.1) 6020 FORMAT(10X,15HGRID DIMENSION=,I3,4H BY,I3//) 6030 FORMAT(10X,23HWINDOW SURROUNDING NO.=,I3,16H MAINT.FACTOR =,F4.2/ 1 10X,12HTAU-NORMAL =,F4.2,20H WINDOW THICKNESS =,F6.2//) 6031 FORMAT(15X,29HWINDOW HAS AN OVERHANG ON TOP/ 1 20X,6HWIDTH=,F4.1,18H DIST.FROM GLASS=,F4.1,15H REFLECTAN 2CE =,F4.2) 6032 FORMAT(15X,31HWINDOW HAS OVERHANGS BOTH SIDES/ 1 20X,6HWIDTH=,F4.1,18H DIST.FROM GLASS=,F4.1,15H REFLECTAN 2CE =,F4.2) 6033 FORMAT(15X,34HWINDOW HAS OVERHANG ON BOTTOM ONLY/ 1 20X,6HWIDTH=,F4.1,18H DIST.FROM GLASS=,F4.1,15H REFLECTAN 2CE =,F4.2) 6034 FORMAT(15X,37HWINDOW HAS OVERHANGS ON TOP AND SIDES/ 1 20X,6HWIDTH=,F4.1,18H DIST.FROM GLASS=,F4.1,15H REFLECTAN 2CE =,F4.2) 6035 FORMAT(15X,38HWINDOW HAS OVERHANGS ON ALL FOUR SIDES/ 1 20X,6HWIDTH=,F4.1,18H DIST.FROM GLASS=,F4.1,15H REFLECTAN 2CE =,F4.2) 6064 FORMAT(10X,47HWINDOW HAS SHEER CURTAIN, SHEER CURTAIN FACTOR=, 1 F6.2//) 6067 FORMAT(10X,40HDIRECT-SUN TRANSMISSION FRACTION, ALPHA ) 6068 FORMAT(15X,6HWINDOW,I3,8H ALPHA=,F5.3/) 6036 FORMAT(15X,18HOVERHANG IS OPAQUE//) 6037 FORMAT(15X,45HOVERHANG IS SEMI-TRANSPARENT, TRANSMISSIVITY=,F4.2 1//) 6038 FORMAT(15X,40HOVERHANG IS TRANSLUCENT, TRANSMISSIVITY=,F4.2//) 6040 FORMAT(5X,15HSURROUNDING NO.,I3/5X,18H******************// 1 10X,23HNUMBER OF OPAQUE WALLS ,I3/) 6041 FORMAT(//2X,37HINSIDE-OVERHANG LUMINANCES FOR WINDOW,I2, 1 15H (FOOT-LAMBERT)/ 2 2X,39H***************************************, 3 15H***************/) 6042 FORMAT(10X,10HTOP ,F7.1) 6043 FORMAT(10X,10HLEFT SIDE ,F7.1/ 1 10X,10HRIGHT SIDE,F7.1) 6044 FORMAT(10X,10HBOTTOM ,F7.1) 6045 FORMAT(10X,10HTOP ,F7.1/ 1 10X,10HLEFT SIDE ,F7.1/ 2 10X,10HRIGHT SIDE,F7.1) 6046 FORMAT(10X,10HTOP ,F7.1/ 1 10X,10HLEFT SIDE ,F7.1/ 2 10X,10HRIGHT SIDE,F7.1/ 3 10X,10HBOTTOM ,F7.1) 6050 FORMAT(1H1////5X,30HDIRECT HORIZONTAL INSOLATION ,F6.2,1X, 1 16HWATTS PER SQ.FT./5X,18HLUMINOUS EFFICACY ,F6.1,1X, 2 15HLUMENS PER WATT/5X,30HDIFFUSE HORIZONTAL INSOLATION , 3 F6.2,1X,16HWATTS PER SQ.FT./5X,18HLUMINOUS EFFICACY , 4 F6.1,1X,15HLUMENS PER WATT/ 5 5X,22HSKY MODEL ,I3/ 6 5X,20HGROUND REFLECTANCE ,F4.2) 6070 FORMAT(//5X,15HWORKING SURFACE,I3/5X,18H******************// 1 10X,4HX0 =,F6.2,6H Y0 =,F6.2,6H Z0 =,F6.2/ 2 10X,4HX1 =,F6.2,6H X2 =,F6.2,6H X3 =,F6.2,5H Y =,F6.2/ 3 10X,7HBETAX =,F6.1,8H PSIX =,F6.1,9H BETAY =,F6.1, 4 8H PSIY =,F6.1/ 5 10X,15HGRID DIMENSION=,I3,4H BY,I3//) 6100 FORMAT(///10X,46HSKY LUMINANCE DETERMINED BY GEOGRAPHICAL DATA / A 15X,24HLATITUDE ,F6.1,5H DEG./ B 15X,24HLONGITUDE ,F6.1,5H DEG./ C 15X,24HTIME ZONE ,I4/ D 15X,24HALTITUDE ,F6.1,7H METERS/ E 15X,24HSKY MODEL ,I4/ F 15X,24HGROUND REFLECTANCE ,F7.2) 6110 FORMAT( 15X,24HFIRST MONTH ,I4/ A 15X,24HLAST MONTH ,I4/ B 15X,24HINCREMENT BETWEEN MONTHS,I4/ C 15X,24HDAY OF MONTH ,I4/ D 15X,24HFIRST TIME OF DAY ,F6.1,5H HRS./ E 15X,24HLAST TIME OF DAY ,F6.1,5H HRS./ F 15X,24HTIME INCREMENT ,F6.1,5H HRS.) 6105 FORMAT(///10X,46HSKY LUMINANCE DETERMINED BY SUN POSITON DATA / A 15X,24HALTITUDE ANGLE ,F6.1,5H DEG./ B 15X,24HAZIMUTH ANGLE ,F6.1,5H DEG./ D 15X,24HALTITUDE OF STATION ,F6.1,7H METERS/ E 15X,24HSKY MODEL ,I4/ F 15X,24HGROUND REFLECTANCE ,F7.2) 6120 FORMAT(15X,45HMONTHLY VALUES FOR TURBIDITY AND COND.WATER../ A 20X,20HMONTH C.W. TURB.) 6121 FORMAT(21X,I3,2F8.4) 6130 FORMAT(//,5X,7HMONTH =,I3,10X,6HTIME =,F5.1,5H HRS./ 1 5X,36H====================================) 6140 FORMAT(5X,20HSUN ALTITUDE ANGLE =,F6.1,/ 15X,19HSUN AZIMUTH ANGLE =,F6.1,1X,27HDEG. OFF SOUTH TOWARDS EAST) 6200 FORMAT(7X,15F6.0) 6250 FORMAT(//10X,35HILLUMINANCE ON A HORIZONTAL SURFACE/ 1 15X,16HSOLAR COMPONENT=,F8.0,12H FOOT-CANDLE/ 2 15X,16HSKY COMPONENT=,F8.0,12H FOOT-CANDLE/ 3 10X,17HZENITH LUMINANCE=,F8.0,1X,14H(FT.-LAMBERTS)) 6300 FORMAT(1H1/15X,47H***********************************************/ * 15X,47H* */ 1 15X,47H* DAYLIGHT ILLUMINATION IN A RECTANGULAR ROOM */ * 15X,47H* */ 2 15X,47H***********************************************// 3 10X,33HROOM HAS A WIDTH (WINDOW WALL) OF,F5.1,6H UNITS/ 4 10X,33H A DEPTH OF,F5.1,6H UNITS/ 5 10X,33H AND A HEIGHT OF,F5.1,6H UNITS/ 6 10X,47HROOM ORIENTATION (FROM OUTWARD WINDOW NORMAL) =, 7 F5.1,15H DEG. OFF SOUTH/) 6302 FORMAT(10X,8HROOM HAS,I2,28H WINDOW(S) IN THE FRONT WALL) 6303 FORMAT(10X,8HROOM HAS,I2,37H WINDOW(S) IN THE CEILING (SKYLIGHTS)) 6310 FORMAT(/10X,10HWINDOW NO.,I2,3H IS,F6.2,14H UNITS WIDE BY,F6.2, 1 11H UNITS HIGH/15X,16HWINDOW CENTER IS,F6.2, 2 30H UNITS TO RIGHT OF WALL CENTER/ 3 31X,F6.2,17H UNITS FROM FLOOR/ 4 15X,24HWINDOW DISCRETIZATION IS,I3,3H BY,I3) 6311 FORMAT(15X,16HWINDOW TYPE IS ",I1,16H" (CLEAR WINDOW)) 6312 FORMAT(15X,16HWINDOW TYPE IS ",I1,34H" (SHEER CURTAIN, CLEAR FRACT 1ION =,F4.2,1H)) 6313 FORMAT(15X,16HWINDOW TYPE IS ",I1,20H" (DIFFUSING WINDOW)) 6320 FORMAT(15X,23HMAINTENANCE FACTOR =,F4.2/ 1 15X,23HWINDOW THICKNESS =,F4.2/ 2 15X,23HWINDOW TRANSMISSIVITY =,F4.2/ 3 15X,23HWINDOW REFLECTANCE =,F4.2) 6331 FORMAT(10X,32HFRONT WALL.. DISCRETIZATION,I3,3H BY,I3, 1 16H, REFLECTANCE =,F4.2) 6332 FORMAT(10X,32HLEFT SIDE.. DISCRETIZATION,I3,3H BY,I3, 1 16H, REFLECTANCE =,F4.2) 6333 FORMAT(10X,32HREAR WALL.. DISCRETIZATION,I3,3H BY,I3, 1 16H, REFLECTANCE =,F4.2) 6334 FORMAT(10X,32HRIGHT SIDE.. DISCRETIZATION,I3,3H BY,I3, 1 16H, REFLECTANCE =,F4.2) 6335 FORMAT(10X,32HCEILING.. DISCRETIZATION,I3,3H BY,I3, 1 16H, REFLECTANCE =,F4.2) 6336 FORMAT(10X,32HFLOOR.. DISCRETIZATION,I3,3H BY,I3, 1 16H, REFLECTANCE =,F4.2) 6337 FORMAT(10X,32HWORKING SURFACE.. DISCRETIZATION,I3,3H BY,I3, 1 13H, ELEVATION =,F5.2,6H UNITS) 6340 FORMAT(/10X,18HSUBJECT BUILDING../ 1 15X,32HBUILDING WIDTH (FRONT WALL) =,F6.1,6H UNITS/ 2 15X,32HBUILDING HEIGHT =,F6.1,6H UNITS/ 3 15X,32HOFFSET TO RIGHT OF ROOM CENTER =,F6.1,6H UNITS/ 4 15X,32H REFLECTANCE =,F7.2) 6350 FORMAT(/10X,13HOBSTRUCTION../ 1 15X,32HOBSTRUCTION WIDTH =,F6.1,6H UNITS/ 2 15X,32HOBSTRUCTION HEIGHT =,F6.1,6H UNITS/ 3 15X,32HOFFSET TO RIGHT OF ROOM CENTER =,F6.1,6H UNITS/ 4 15X,32HDISTANCE FROM SUBJECT BUILDING =,F6.1,6H UNITS/ 5 15X,32H OBSTRUCTION REFLECTANCE =,F7.2) 6560 FORMAT(//2X,69HLIGHT EXCHANGE FACTORS FROM SURFACE I TO SURFACE J 1IN SURROUNDING NO.,I2,/2X,71H************************************* 2**********************************/) 6570 FORMAT(/7X,1HJ,15I8) 6572 FORMAT(4X,1HI) 6580 FORMAT(2X,I3,5X,15F8.4/) 6710 FORMAT(1H1//3X,71H************************************************ 1***********************/ 2 3X,71H* DATA FOR SURFACE NODES; I=SURFACE, K=NODE, X,Y, 3Z=COORDINATES (FT) */ 4 3X,71H* A=AREA (SQFT), L=LUMINANC 5E (FT-LAMBERT) */ 6 3X,71H*************************************************** 7********************) 6715 FORMAT(1H1//3X,71H************************************************ 1***********************/ 2 3X,71H* DATA FOR WORKING SURFACE NODES; I=SURFACE, K=NODE 3-NO. */ 4 3X,71H* X,Y,Z=COORDINATES, 5 */ ) 6716 FORMAT( 3X,71H* S=ILLUMINANCE FROM EXTERNAL SOURCES (FOOT-CA 7NDLE) */ 8 3X,71H* R=INTERNAL REFLECTED COMPONENT (FOOT-CANDLE) 9 */ A 3X,71H* I=TOTAL ILLUMINANCE (FOOT-CANDLE), D=DAYLIGH BT FACTOR (PERCENT) */ 6 3X,71H*************************************************** 7********************) 6717 FORMAT(//10X,53HX = DISTANCE FROM FRONT WALL (PARALLEL TO FRONT WA 1LL)/ 2 10X,64HY = DISTANCE FROM ROOM CENTER LINE (PERPENDICULAR 3TO FRONT WALL)/ 4 14X,37HPOSITIVE VALUES (LEFT OF CENTER LINE)/ 5 14X,38HNEGATIVE VALUES (RIGHT OF CENTER LINE)/ 6 10X,24HZ = HEIGHT OF WORK PLANE) 6720 FORMAT(//3X,2H*K,12I10) 6730 FORMAT(3X,2HI*/I4,3H X,12F10.2) 6740 FORMAT(6X,1HY,12F10.2) 6750 FORMAT(6X,1HZ,12F10.2) 6760 FORMAT(6X,1HA,12F10.2) 6770 FORMAT(6X,1HL,12F10.2) 6773 FORMAT(6X,1HS,12F10.2) 6774 FORMAT(6X,1HR,12F10.2) 6775 FORMAT(6X,1HI,12F10.2) 6780 FORMAT(6X,1HD,12F10.2) 7010 FORMAT(//5X,41HLUMINANCE OF SURFACE I IN SURROUNDING NO.,I2,1X,14H 1(FT.-LAMBERTS)/5X,58H============================================= 2=============) 7020 FORMAT(/5X,1HI,15I6) 8000 FORMAT(2X,33HTHERE IS TOO MUCH DIRECT SUNLIGHT/ 1 2X,35HAS COMPARED TO TOTAL HEMISPHERICAL,/ 2 2X,20HEXECUTION TERMINATED) 8002 FORMAT(32HNO. OF WINDOWS EXCEEDS MAXIMUM 5) 8005 FORMAT(33HILLEGAL REFLECTANCE FOR SURFACE ,I2, 2 14H IN ENCLOSURE ,I2) 8007 FORMAT(38HNO. OF WORK SURFACES EXCEEDS MAXIMUM 3) 8008 FORMAT(35HNO. OF ENCLOSURES EXCEEDS MAXIMUM 2) 8010 FORMAT(41HNO. OF INDOOR SURFACES EXCEEDS MAXIMUM 30) 8011 FORMAT(44HNO. OF NODES EXCEEDS 400 FOR INDOOR SURFACE ,I2) 8012 FORMAT(40HILLEGAL REFLECTANCE FOR INDOOR SURFACE ,I2) 8014 FORMAT(44HILLEGAL SURFACE OCCUPYING CUTOUT IN SURFACE ,I2) 8020 FORMAT(29HILLEGAL ENCLOSURE FOR WINDOW ,I2) 8022 FORMAT(34HILLEGAL TRANSMISSIVITY FOR WINDOW ,I2) 8024 FORMAT(38HILLEGAL MAINTENANCE FACTOR FOR WINDOW ,I2) 8026 FORMAT(27HILLEGAL GROUND REFLECTANCE ) 8030 FORMAT(54HTHICKNESS OF CONDENSABLE WATER OUT OF RANGE FOR MONTH , 1 I2) 8032 FORMAT(45HTURBIDITY COEFFICIENT OUT OF RANGE FOR MONTH ,I2) 8034 FORMAT(41HILLEGAL OVERHANG REFLECTANCE FOR WINDOW ,I2) 8036 FORMAT(43HILLEGAL OVERHANG TRANSMISSIVITY FOR WINDOW ,I2) 8100 FORMAT(6F10.3) 8200 FORMAT(20I4) END FUNCTION C(K,NK) C THIS FUNCTION CALCULATES NEWTON-COTES COEFFICIENTS FOR NUMERICAL C QUADRATURE C=1. IF(NK.EQ.1)RETURN I=(K+1)/2-K/2 A=1+2*(1-I)+((K+100)/102)-(K/NK) C=A/(3.*(NK-1.)) RETURN END SUBROUTINE SECT(X,Y,K,J,KUT) C THIS ROUTINE CHECKS WHETHER A RAY GOING FROM X INTO DIRECTION Y CUTS C SURFACE J IN ENCLOSURE K FROM ITS FRONT SIDE (KUT=1) OR NOT (KUT=0) REAL X(3),Y(3),X0(4,30,7),LL(4,30,9) COMMON /BLK1/ X0,LL,XDP,YDP C C INTERSECTION POINT IN LOCAL COORDINATES BN=0. BD=0. DO 1 I=1,3 BN=BN+(X(I)-X0(K,J,I))*LL(K,J,I+6) 1 BD=BD+ Y(I)*LL(K,J,I+6) BD=BD+1E-10+BD*1E-5 B=BN/BD XDP=0. YDP=0. DO 2 I=1,3 CC=X(I)-X0(K,J,I)-B*Y(I) XDP=XDP+CC*LL(K,J,I) 2 YDP=YDP+CC*LL(K,J,I+3) KUT=0 C WITHIN CONFINES OF SURFACE? ETA=YDP/X0(K,J,7) IF(ETA.LT.0..OR.ETA.GT.1.) RETURN XI=(XDP-ETA*X0(K,J,5))/(X0(K,J,4)-(X0(K,J,4)-X0(K,J,6)+ 1 X0(K,J,5))*ETA) IF(XI.LT.-1.E-5.OR.XI.GT.1.00001) RETURN C IS IT HITTING TOP OF SURFACE? KUT=1-(Y(1)*LL(K,J,7)+Y(2)*LL(K,J,8)+Y(3)*LL(K,J,9))*1E-4 RETURN END SUBROUTINE CONFAC(J,NT,F) C THIS ROUTINE CALCULATES EXCHANGE FACTORS FOR ENCLOSURE J WITH NT C SURFACES, INCLUDING EXCHANGE FACTORS TO SKY WITH LUMINANCE VARIATION REAL X0(4,30,7),LL(4,30,9),F(4,30,30),X(3),Y(3), 1 A(4,30),RS(4),LPZ INTEGER NW(6),IOB(4,30,30),KJOB(4,30,30) INTEGER*4 ISEED COMMON /BLK1/ X0,LL,XDP,YDP COMMON /BLK2/ NW,PI,IOB,ICLD,ISEED,NBUNDL,A,KJOB NT1=NT+1 DO 10 I=1,NT DO 5 IR=1,NT1 5 F(J,I,IR)=0. XY=(X0(J,I,6)-X0(J,I,5))/X0(J,I,4) ID1=1 IF(ABS(XY-1).LT.1E-3)ID1=0 AX=ID1*XY C CALCULATION OF COORDINATES AND DIRECTION COSINES OF LIGHT BUNDLE DO 60 N=1,NBUNDL TT=1. C ************************ RANDOM NUMBER GENERATOR ********************** DO 115,INDEX=1,4 115 RS(INDEX)=RANDOM(ISEED) YY=(1-SQRT(1-(1-XY*XY)*RS(1)))/(1-AX)+(1-ID1)*RS(1) XX=X0(J,I,5)*YY+X0(J,I,4)*(1-(1-XY)*YY)*RS(2) SINB=SQRT(RS(3)) COSB=SQRT(1.-RS(3)) THET=2.*PI*RS(4) COST=COS(THET) SINT=SIN(THET) DO 30 L=1,3 X(L)=X0(J,I,L)+LL(J,I,L)*XX+LL(J,I,L+3)*YY*X0(J,I,7) Y(L)=(LL(J,I,L)*COST+LL(J,I,L+3)*SINT)*SINB+LL(J,I,L+6)*COSB 30 CONTINUE C C DETERMINATION OF SURFACE INTERSECTED BY LIGHT BUNDLE C DO 80 IR=1,NT IF(IOB(J,I,IR)) 80,70,50 50 IMARK=KJOB(J,I,IR) CALL SECT(X,Y,J,IMARK,KUT) IF(KUT.EQ.1) GOTO 60 70 IMARK=IR CALL SECT(X,Y,J,IR,KUT) IMARK=IR IF(KUT.EQ.1) GOTO 60 80 CONTINUE C C EVALUATION OF F I-J IF BUNDLE INTERSECTS SKY, (OR GROUND IF Y POINTS BELOW HORIZ) C IMARK=NT1 TT=RLPZ(ICLD,Y) 60 F(J,I,IMARK)=F(J,I,IMARK)+TT DO 130 IR=1,NT1 130 F(J,I,IR)=F(J,I,IR)/NBUNDL 10 CONTINUE C C AVERAGING OF F I-J BASED UPON SURFACE AREA C DO 220 I=1,NT DO 220 IR=I,NT WF=1-EXP(-A(J,I)/A(J,IR))+EXP(-A(J,IR)/A(J,I)) AF=((2-WF)*A(J,I)*F(J,I,IR)+WF*A(J,IR)*F(J,IR,I))/2. F(J,I,IR)=AF/A(J,I) 220 F(J,IR,I)=AF/A(J,IR) RETURN END FUNCTION TAUB(I,CB) C THIS ROUTINE CALCULATES THE TRANSMISSIVITY OF WINDOW I FOR DIRECTION C CB (COSINE OF POLAR ANGLE) REAL AK2(10),TAU(10) COMMON /BLK3/ AK2,TAU,NC TAUB=1.018*AK2(I)*TAU(I)*CB*(1.+(1.-CB*CB)**1.5) IF (I.GT.NC) TAUB=TAU(I) RETURN END FUNCTION RLPZ(ICLD,Y) C THIS ROUTINE CALCULATES SKY-LUMINANCE (NORMALIZED BY ZENITH LUMINANCE) C FOR ANY GIVEN DIRECTION AND SKY CONDITION REAL Y(3),LS(3) COMMON /BLK4/ BETAS,LS,AUX(3),GRNDL IF (Y(3).GT..001) GOTO 10 C HITS GROUND RLPZ=GRNDL RETURN 10 GOTO(1,2,3,2),ICLD C OVERCAST SKY 1 IY=.999+Y(3) RLPZ=AUX(1)+IY*Y(3)*AUX(2) RETURN C CLEAR SKY 2 CD=Y(1)*LS(1)+Y(2)*LS(2)+Y(3)*LS(3) RLPZ=(1.-EXP(-.32/Y(3)))*(.91+10*EXP(-3*ACOS(CD))+.45*CD*CD) 1 /AUX(3) RETURN C UNIFORM SKY 3 RLPZ=1. RETURN END SUBROUTINE MATINV(A,N,B,M,DETERM) DIMENSION IPIVOT(30),A(30,30),B(30,1),INDEX(30,2),PIVOT(30) EQUIVALENCE (IROW,JROW), (ICOLUM,JCOLUM), (AMAX, T, SWAP) DETERM=1.0 DO 20 J=1,N 20 IPIVOT(J)=0 DO 550 I=1,N AMAX=0.0 DO 105 J=1,N IF (IPIVOT(J)-1) 60, 105, 60 60 DO 100 K=1,N IF (IPIVOT(K)-1) 80, 100, 740 80 IF (ABS(AMAX)-ABS(A(J,K))) 85, 100, 100 85 IROW=J ICOLUM=K AMAX=A(J,K) 100 CONTINUE 105 CONTINUE IF(AMAX) 110,800,110 110 IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1 IF (IROW-ICOLUM) 140, 260, 140 140 DETERM=-DETERM DO 200 L=1,N SWAP=A(IROW,L) A(IROW,L)=A(ICOLUM,L) 200 A(ICOLUM,L)=SWAP IF(M) 260, 260, 210 210 DO 250 L=1, M SWAP=B(IROW,L) B(IROW,L)=B(ICOLUM,L) 250 B(ICOLUM,L)=SWAP 260 INDEX(I,1)=IROW INDEX(I,2)=ICOLUM PIVOT(I)=A(ICOLUM,ICOLUM) DETERM=DETERM*PIVOT(I) A(ICOLUM,ICOLUM)=1.0 DO 350 L=1,N 350 A(ICOLUM,L)=A(ICOLUM,L)/PIVOT(I) IF(M) 380, 380, 360 360 DO 370 L=1,M 370 B(ICOLUM,L)=B(ICOLUM,L)/PIVOT(I) 380 DO 550 L1=1,N IF(L1-ICOLUM) 400, 550, 400 400 T=A(L1,ICOLUM) A(L1,ICOLUM)=0.0 DO 450 L=1,N 450 A(L1,L)=A(L1,L)-A(ICOLUM,L)*T IF(M) 550, 550, 460 460 DO 500 L=1,M 500 B(L1,L)=B(L1,L)-B(ICOLUM,L)*T 550 CONTINUE DO 710 I=1,N L=N+1-I IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630 630 JROW=INDEX(L,1) JCOLUM=INDEX(L,2) DO 705 K=1,N SWAP=A(K,JROW) A(K,JROW)=A(K,JCOLUM) A(K,JCOLUM)=SWAP 705 CONTINUE 710 CONTINUE 740 RETURN 800 DETERM = 0. RETURN END SUBROUTINE ZENLCD(P,W,BETA,H,ZENL) C THIS SUBROUTINE CALCULATES ZENITH LUMINANCE BASED ON C LIEBELT'S EQUATION WHEN SOLAR ALTITUDE IS LESS THAN C 60. WHEN SOLAR ALTITUDE IS GREATER THAN 60, THE ZENITH C LUMINANCE IS EXTRAPOLATED TO RESULT INTEGRATED HORIZONTAL C ILLUMINANCE FOLLOWS SINE CURVE. 4/10/83 JJK. PI=4.*ATAN(1.) HD=H*180./PI T=TUR(HD,W,BETA) ZENL=((1.34*T-3.46)*TAN(H)+0.1*T+0.9)*291.87 IF(HD.LT.60.)RETURN ZEN60=((1.34*T-3.46)*TAN(1.0472)+0.1*T+0.9)*291.87 P60=SKYSUM(0.5236) ZENL=1.1547*ZEN60*P60*SIN(H)/P RETURN END SUBROUTINE SOLAR (DLAT,LNG,TZN,M,TIME,ISUN,THETA,PHI) C TO FIND THE POSITION OF THE SUN C INPUT C DLAT - LATITUDE C LNG - LONGITUDE C TZN - TIME ZONE(HOURS BEHIND GREENWICH MEAN TIME C 4- ATLANTIC C 5- EASTERN C 6- CENTRAL C 7- MOUNTAIN C 8- PACIFIC C M - MONTH(01-12) C TIME - HOURS FROM MIDNIGHT C C OUTPUT C THE ALTITUDE AND AZIMUTH ANGLES OF THE SUN. REAL LAT, LNG INTEGER TZN DIMENSION DEC(12), ET(12) C DATA (DEC(J), J=1,12)/ -0.3491, -0.1885, 0.0 , 0.2025, $ 0.3491, 0.4093, 0.3595, 0.2147, $ 0.0 , -0.1833, -0.3456, -0.4093/ C DATA (ET(J), J=1,12)/ -0.190, -0.230, -0.123, 0.020, $ 0.060, -0.025, -0.103,-0.051, $ 0.113, 0.255, 0.235, 0.033/ DATA PI/3.1415926/ DATA CNVRT/0.017453295/ C WHEN SUN IS DOWN THERE IS NO RADIATION LAT = DLAT*CNVRT SPAN = ACOS(-TAN(LAT)*TAN(DEC(M))) HANG =-(15.0*(TIME - 12.0 + TZN + ET(M)) - LNG)*CNVRT CHANG = COS(HANG) IF(ABS(HANG) .GT. ABS(SPAN)) THEN ISUN = 0 C WHEN SUN IS UP, FIND ITS ORIENTATION ELSE ISUN = 1 COSZ = SIN(LAT)*SIN(DEC(M)) + $ COS(LAT)*COS(DEC(M))*COS(HANG) COSW = COS(DEC(M))*SIN(HANG) TEST = 1.0 - COSZ**2 + COSW**2 ARG = TAN(DEC(M))/TAN(LAT) IF (TEST .LE. 0.0) THEN COSS = 0.0 ELSEIF (CHANG .LE. ARG) THEN COSS = -SQRT(TEST) ELSE COSS = SQRT(TEST) ENDIF THETA = ASIN(COSZ) CHECK = COSW/COS(THETA) IF (CHECK .GE. 1.0) THEN ANG = PI/2.0 ELSEIF (CHECK .LE. -1.0) THEN ANG = -PI/2.0 ELSE ANG = ASIN(CHECK) ENDIF IF (COSS .GT. 0.0) THEN PHI = ANG ELSE PHI = PI - ANG ENDIF ENDIF RETURN END FUNCTION E(DATE,IWMOD) REAL OMEGA DATA OMEGA/.01715/ A = OMEGA*DATE IF(IWMOD.NE.3) GOTO 100 E=128.8 RETURN 100 E = 128.82 + 4.248*COS(A) + 0.0825*COS(2.0*A) $-0.00043*COS(3.0*A) + 0.169*SIN(A) + 0.00914*SIN(2.0*A) $+0.01726*SIN(3.0*A) RETURN END FUNCTION ABAR(T,BETA) DIMENSION A(3),B(3) DATA (A(I),I=1,3)/0.1512,0.1656,0.2021/ DATA (B(I),I=1,3)/0.0262,0.0215,0.0193/ INDEX=3 IF (BETA .LE. .15) INDEX=2 IF (BETA .LE. .075) INDEX=1 ABAR = A(INDEX) - T*B(INDEX) RETURN END FUNCTION TUR(HD,W,BETA) TUR=(HD + 85.0)/(39.5*EXP(-W) + 47.4) + 0.1 $+ (16.0+0.22*W)*BETA RETURN END FUNCTION AIR(H,ALT) AIR = SIN(H) + 0.15*(H*57.29578+3.885)**(-1.253) AIR = 1/AIR AIR = AIR*(1.0-1.E-4*ALT) RETURN END FUNCTION DIR(MONTH,IDAY,BETA,W,ALT,H) DATA CFAC/92.904/ C CFAC CONVERTS FROM KLX TO FT-CANDLES HD = H*57.2958 T = TUR(HD,W,BETA) A = ABAR(T,BETA) AMAS = AIR(H,ALT) DATE = FLOAT(JDATE(MONTH,IDAY)) DIR = CFAC*E(DATE,IWMOD)*EXP(-A*AMAS*T)*SIN(H) RETURN END FUNCTION JDATE(MONTH,IDAY) DIMENSION NDAY(12) DATA (NDAY(I),I=1,12)/31,28,31,30,31, $30,31,31,30,31,30,31/ JDATE = 0 IF (MONTH .EQ. 1) GO TO 200 M1=MONTH-1 DO 100 I=1,M1 JDATE = JDATE + NDAY(I) 100 CONTINUE 200 JDATE = JDATE + IDAY RETURN END FUNCTION SKYSUM(H) C -- THIS SUBROUTINE CALCULATES SKY LUMINANCE INTEGRATION FOR C-- CLEAR SKY PI=4.*ATAN(1.) P1=.27385*(.91+10.*EXP(-3*H)+.45*COS(H)**2) P=0. CBS=COS(H) SBS=SIN(H) DO 100 K=1,21 ETAK=(K-1.)/20. RT=SQRT(1.-ETAK**2) P2=0. IF(K.GT.1)P2=(1.-EXP(-.32/ETAK))*C(K,21)*ETAK DO 100 L=1,21 ETA=(L-1.)/20. CDEL=CBS*ETAK+SBS*RT*COS(PI*ETA) 100 P=P+P2*(.91+10.*EXP(-3*ACOS(CDEL))+.45*CDEL**2)*C(L,21) SKYSUM=2./P1*P RETURN END C FUNCTION RANDOM(ISEED) INTEGER*4 ISEED ISEED=DMOD(16807.0D0 * ISEED,2147483647.0D0) RANDOM=ISEED * 4.6566128752458D-10 RETURN END SUBROUTINE LOGO(N) C THIS ROUTINE PRINTS THE SUPERLITE NOTICE TO THE UNIT C SPECIFIED AS 'N' C CURRENTLY N=6 SENDS NOTICE TO SCREEN ON VAX. FOR BATCH C MACHINES THIS CAN BE SET TO THE APPROPRIATE DEVICE WRITE(*,100) 100 FORMAT(///////// *6X,' S U P E R L I T E 1 . 0'/ *// *6X,' DEVELOPED BY:'/ */ *6X,' WINDOWS AND DAYLIGHTING GROUP'/ *6X,' LAWRENCE BERKELEY LABORATORY'/ *6X,' UNIVERSITY OF CALIFORNIA'/ *6X,' BERKELEY, CA 94720'/ */// *6X,' COPYRIGHT 1985 REGENTS OF THE UNIVERSITY OF CALIFORNIA' *//) END