PROGRAM MAIN 这是单目标遗传算法.doc

PROGRAM MAIN !这是单目标遗传算法INTEGER POPSIZE !种群规模INTEGER NUM(100)!染色体选择时染色体编号记录数组INTEGER NUMBER !记录选择时染色体编号INTEGER I !循环变量INTEGER J !循环变量INTEGER K !循环变量INTEGER L !循环变量INTEGER N !染色体内基因数（目标函数自变量数）INTEGER YES(100) !约束函数满足与否标志（0为不满足，1为满足）INTEGER GEN !遗传代数!INTEGER IDUM !随机数种子COMMON/X0/X(100)!基因数组COMMON/Y0/Y(100)!遗传操作时基因交叉、变异用数组REAL V(1000,100)!初始染色体数组(1000代表染色体数量，100代表染色体内基因数量)COMMON/RANSETI/ V1(1000,100)! 初始染色体排序后数组COMMON/RANSETI/ V2(1000,100)!旋转轮赌排序后染色体数组COMMON/RANSETI/ V3(1000,100)!交叉操作后染色体数组COMMON/RANSETI/ V4(1000,100)!变异操作后染色体数组COMMON/R/R(10000) !R(I)COMMON/P/P(10000) !P(I)COMMON/PSAI/PSAI(10000) !(I)COMMON/OMG/OMG(10000) !外荷载与竖直线的夹角COMMON/U/U(100000)COMMON/Q0/Q0(10000) !外荷载COMMON/N1/N1 !土条划分数目COMMON/F0/F0 !初始安全系数COMMON/LMD0/LMD0 !初始COMMON/LMD/LMD !计算的COMMON/FS/FS !计算的安全系数COMMON/SP/SPCOMMON/SP1/SP1COMMON/SP2/SP2COMMON/MP/MPCOMMON/M1/M1COMMON/M2/M2COMMON/XDING/XDING,YDING !读入坡顶的坐标值COMMON/N3/N3 !底部划分条块数COMMON/N4/N4 !中部划分条块数COMMON/N5/N5 !顶部划分条块数COMMON/N6/N6 !第2种滑面顶部划分条块数1COMMON/N7/N7 !第2种滑面顶部划分条块数2COMMON/IDUM/IDUM !随机数种子COMMON/C1/C1(1000)COMMON/FAI1/FAI1(1000)COMMON/GAMA1/GAMA1(1000)COMMON/NB/NB !土坡表面几何特征点数NBCOMMON/NS/NS !土层数NSCOMMON/EQH/EQH !水平地震系数EQHCOMMON/XU/XU(10000)!地表特征点水平坐标COMMON/YU/YU(1000,10000) !地表特征点下各土层界面的垂直坐标COMMON/WL/WL(10000) !地表特征点下地下水或浸润线的垂直坐标COMMON/QA/QA(10000) !地表特征点间的荷载1COMMON/ALF/ALF(10000) !各条块的倾角COMMON/W/W(10000) !各条块的重度COMMON/C/C(10000) !各条块的粘聚力COMMON/FAI/FAI(10000)!各条块的内摩擦角COMMON/U/U(10000)!水压力COMMON/WQH/WQH(10000)!水平地震力COMMON/WQV/WQV(10000)!竖直地震力COMMON/X1/X1(10000)!滑动面控制点横坐标COMMON/Y1/Y1(10000)!滑动面控制点纵坐标COMMON/E/E(10000) !E(I)COMMON/B/B(10000) !每一条块的宽度COMMON/H1/H1(10000) !条块上下中心点的距离COMMON/XCENTER/XCENTER(10000)COMMON/YCENTER/YCENTER(10000)COMMON/FAII/FAII(10000)COMMON/PS/PS(10000)COMMON/DILMD/DILMD !底部方向向量值!COMMON/DINGLMD/DINGLMD !顶部方向向量值COMMON/NTIAO/NTIAO !判断属于哪种分条模式COMMON/IYES/IYES !分条模式成功与否判断COMMON/LMDL/LMDL(100) !储存随机产生的基因数组COMMON/LEFT2/LEFT2(100) !基因下限COMMON/RIGHT2/RIGHT2(100) !基因上限COMMON/LMDBEST/LMDBEST(100)COMMON/XIANZHI1/XIANZHI1 !构建滑裂面的限值1COMMON/XIANZHI2/XIANZHI2 !构建滑裂面的限值2COMMON/XIANZHI3/XIANZHI3 !构建滑裂面的限值3COMMON/XIANZHI4/XIANZHI4 !构建滑裂面的限值4COMMON/XIANZHI5/XIANZHI5 !安全系数限制COMMON/LEFT1/LEFT1(100) !基因下限COMMON/RIGHT1/RIGHT1(100) !基因上限REAL F(100)!目标问题函数值记录数组REAL Q(100)!评价函数值记录数组REAL FATHER(1000,100)!遗传操作时父代存储数组REAL CHILD(1000,100)! 遗传操作时子代存储数组REAL TL(100),YLR(100),BLR,RLR !非均匀变异的参数T,Y,D,BREAL VLEFT(10000,50) !记录染色体的下限REAL VRIGHT(1000,50) !记录染色体的上限REAL FF(10000) !,VC(10000,50)!记录随机搜索的安全系数及染色体INTEGER POPSIZE2,NUM2(10000),YES2(10000) !初始随机搜索次数REAL VV(10000,50),VV1(10000,50)!记录随机搜索的染色体REAL MIN !过渡量REAL TEMP !过渡量REAL BEST !每代染色体中最优目标函数值记录INTEGER BEST1 !每代染色体中最优目标函数值记录REAL BESTG(100) !每代染色体中最优目标函数值对应染色体记录REAL LEFT1(100),RIGHT1(100),LEFT2REAL M !变异时方向系数REAL FAVG1,FAVG2 ! 平均适应值REAL FMAX2 ! 最大适应值REAL FMIN1,FMIN2 ! 最小适应值REAL FIT1(1000),FIT2(1000) !适应值数组REAL TEMP4,TEMP5 !临时变量REAL PC(1000) !交叉数组REAL PM(1000) !变异数组2REALREALREALREALPC1,PC2 !最大及最小交叉率PM1,PM2 !最大及最小变异率M1,M2,MPLMD,LMD0,LMDL,LMDBESTPARAMETER(PI=3.1415926)OPEN(3,FILE=滑面图.SCR)OPEN(4,FILE=输入数据.DAT,STATUS=OLD)OPEN(5,FILE=结果.DAT)OPEN(6,FILE=滑面参数.DAT)OPEN(7,FILE=遗传结果.DAT)OPEN(8,FILE=适应度结果.DAT)READ(4,*)GEN,POPSIZE,N,IDUM,PC1,PC2,PM1,PM2,BLR!读入遗传代数，种群规模，染色体内基因数,随机数种子,交叉操作概率,变异操作概率,变异指数READ(4,*)EPXL1,EPXL2!10!WRITE(*,*)EPXL1,EPXL2READ(4,*)N3,N4,N5,N6,N7 !读入各控制段要划分的条块数目WRITE(*,*)N2,N3,N4,N5,N6,N7READ(4,*)IDUMWRITE(*,*)IDUMREAD(4,*)NB,NS !输入:土坡表面几何特征点数NB,土层数NS,土条数N1READ(4,*)EQHREAD(4,*)(XU(J),J=1,NB) !输入地表特征点水平坐标:,X1,X2,.Xnb=DO 10 I=1,NSWRITE(*,(2X,A,I2,A)输入第,I,层:Y1,Y2.Ynb=READ(4,*)(YU(I,J),J=1,NB)CONTINUEREAD(4,*)(WL(J),J=1,NB)!输入地表特征点下地下水或浸润线的垂直坐标:!,WL1,WL2.WLnb=READ(4,*)(QA(J),J=1,NB) !输入地表特征点间的荷载:,QA1,QA2.QAnb=READ(4,*)(GAMA1(I),I=1,NS)!读入各层的重度READ(4,*)(C1(I),I=1,NS)!读入各层的粘聚力READ(4,*)(FAI1(I),I=1,NS)!读入各层的内摩擦角READ(4,*)LEFTDI,RIGHTDI !读入底部方向向量的上下限WRITE(*,*)LEFTDI,RIGHTDIREAD(4,*)LEFTDING,RIGHTDING !读入顶部方向向量的上下限READ(4,*)NTIAO !读入分条模式形式READ(4,*)XDING,YDING !读入坡顶坐标READ(4,*)SP1,SP2,M1,M2READ(4,*)XIANZHI1,XIANZHI2,XIANZHI3,XIANZHI4READ(4,*)XIANZHI5,POPSIZE2DO I=1,NREAD(4,*)LEFT1(I),RIGHT1(I)ENDDOWRITE(*,*)SP1,SP2,M1,M2CLOSE(4)STOPDO I=1,NSFAI1(I)=FAI1(I)*PI/180.0ENDDODO 30 I=1,POPSIZE2WRITE(*,*)I=,INUM2(I)=I315!20!25!30!X(1)=LEFTDI+RAN3(IDUM)*(RIGHTDI-LEFTDI)X(N)=LEFTDING+RAN3(IDUM)*(RIGHTDING-LEFTDING)WRITE(*,*)X(1),X(N)DILMD=X(1)DINGLMD=X(N)WRITE(*,*)DILMD,DINGLMDCALL RENYIHUAMIANWRITE(*,*)IYES=,IYESIF(IYES.EQ.0)STOPIF(IYES.EQ.0)GOTO 15stopWRITE(8,*)种群序数=,IDO 20 J=1,NX(J)=LMDL(J)VLEFT(I,J)=LEFT2(J)VRIGHT(I,J)=RIGHT2(J)WRITE(*,*)X(J)WRITE(8,*)VLEFT(I,J),VRIGHT(I,J)CONTINUEWRITE(*,*)111111CALL GX1(YES2(I)!调用约束函数WRITE(*,*)111WRITE(*,*)12IF (YES2(I).EQ.1) THENWRITE(*,*)YES(,I,)=,YES(I)WRITE(*,*)X(1),X(2),X(3)CALL MUBIAO(FF(I)!若基因可行则调用目标函数IF (FF(I).GT.50) GOTO 15WRITE(*,*)这里错了？WRITE(*,*)FF(,I,)=,FF(I)WRITE(*,*)(X(J),J=1,N)WRITE(*,*)13DO 25 J=1,NVV(I,J)=X(J)!若基因可行则记录于对应染色体中WRITE(*,*)VV(I,J)CONTINUEELSEGOTO 15END IFSTOPWRITE(5,*)(X(J),J=1,N)WRITE(*,*)F(,I,)=,F(I)CONTINUESTOP!初始染色体组组排序，按从小到大排列DO I=1,POPSIZE2MIN=FF(I)WRITE(*,*)MINNUMBER=NUM2(I)DO 45 J=I,POPSIZE2IF (FF(J).LT.MIN) THENTEMP=MIN4MIN=FF(J)FF(J)=TEMPTEMP=NUMBERNUMBER=NUM2(J)NUM2(J)=TEMPENDIF45 CONTINUEFF(I)=MINNUM2(I)=NUMBERDO 65 J=1,NVV1(I,J)=VV(NUM2(I),J)! WRITE(*,*)VV1(I,J)65 CONTINUE !生成V1(初始染色体排序后数组)! WRITE(5,*)I,FF(I),(VV1(I,J),J=1,N)ENDDO! STOPDO 1001 J=1,NLEFT1(J)=VLEFT(1,J)RIGHT1(J)=VRIGHT(1,J)DO 1002 I=1,POPSIZE2IF(LEFT1(J).GT.VLEFT(I,J)THENLEFT1(J)=VLEFT(I,J)ENDIFIF(RIGHT1(J).LT.VRIGHT(I,J)THENRIGHT1(J)=VRIGHT(I,J)ENDIF1002 CONTINUE1001 CONTINUE!记录最小值DO 1003 I=1,POPSIZEF(I)=FF(I)DO 1004 J=1,N!10041003!56WRITE(*,*)VV1(I,J)V1(I,J)=VV1(I,J)NUM(I)=NUM2(I)WRITE(*,*)NUM(I)WRITE(*,*)V1(I,J)CONTINUECONTINUEBEST=F(1)BEST1=1DO 56 I=1,NBESTG(I)=V1(1,I)WRITE(*,*)BESTG(I)CONTINUE!输出初始最优目标函数值WRITE(5,1000)BEST1000 FORMAT(1X,初始最优目标函数值=,3X,F16.8)WRITE(5,*)(BESTG(I),I=1,N)!1005 FORMAT(1X,初始最优染色体为,3X,F10.6,3X,F10.6,3X,F10.6)WRITE(5,*)RIGHTDI=RIGHT1(1)5!遗传循环DO 10000 L=1 , GENWRITE(*,*)L=,LWRITE(5,1025)L1025 FORMAT(1X,第(,I4,)代数据)!生成评价函数值Q(1)=0TEMP=0DO 75 I=1 , POPSIZECALL EVAL(I,Q(I+1)75807374!7685!88!104090!CONTINUEDO 80 I=1 , POPSIZEQ(I+1)=Q(I+1)/Q(POPSIZE+1)CONTINUE!生成评价函数值结束!适应度函数及其尺度变换TEMP4=0.0DO 73 I=1,POPSIZETEMP4=TEMP4+F(I)CONTINUEFAVG1=TEMP4/POPSIZEDO 74 I=1,POPSIZEFIT1(I)=1.0/(1+F(POPSIZE)+F(I)CONTINUEFMIN1=FIT1(POPSIZE)-0.000001DO 76 I=1,POPSIZEFIT1(I)=FAVG1*(FIT1(I)-FMIN1)/(FAVG1-FMIN1)WRITE(*,*)FIT1(,I,)=,FIT1(I)CONTINUE!适应度函数及尺度变换结束!旋转轮赌选择DO 90 I=1,POPSIZETEMP=RAN3(IDUM)DO 88 J=1,POPSIZEIF (Q(J).LT.TEMP.AND.TEMP.LE.Q(J+1)THENDO 85 K=1 , NV2(I,K)=V1(J,K)FIT2(I)=FIT1(J)CONTINUEWRITE(*,*)FIT2(,I,)=,FIT2(I)ENDIFCONTINUEWRITE(5,1040)I,(V2(I,J),J=1,N)FORMAT(1X,V2(,I2,)=,3X,F8.4,3X,F8.4,3X,F8.4)CONTINUE !旋转轮赌选择结束!计算各染色体的交叉概率TEMP5=0.0FMIN2=FIT2(POPSIZE)FMAX2=FIT2(1)WRITE(8,*)FMIN2=,FMIN2WRITE(8,*)FMAX2=,FMAX2DO 91 I=1,POPSIZE6TEMP5=TEMP5+FIT2(I)IF(FMIN2.GT.FIT2(I)FMIN2=FIT2(I)IF(FMAX2.LT.FIT2(I)FMAX2=FIT2(I)91!9295100!130!105!115CONTINUEWRITE(8,*)FMIN2=,FMIN2WRITE(8,*)FMAX2=,FMAX2FAVG2=TEMP5/POPSIZEDO 92 I=1,POPSIZEIF (FIT2(I).LT.FAVG2.OR.FMAX2.EQ.FAVG2)THENPC(I)=PC1ELSEPC(I)=PC1-(PC1-PC2)*(FIT2(I)-FAVG2)/(FMAX2-FAVG2)ENDIFWRITE(*,*)PC(,I,)=,PC(I)CONTINUE!选择交叉用染色体NUMBER=0DO 100 I=1,POPSIZETEMP=RAN3(IDUM)IF (TEMP.LT.PC(I) THENNUMBER=NUMBER+1NUM(NUMBER)=IDO 95 J=1 , NFATHER(NUMBER,J)=V2(I,J)CONTINUEENDIFCONTINUE !交叉父代选择结束K=INT(NUMBER/2)WRITE(*,*)K=,KDO 110 I= 1 , KWRITE(*,*)I=,ITEMP=RAN3(IDUM)DO 105 J=1,NX(J)=TEMP*FATHER(2*I-1,J)+(1-TEMP)*FATHER(2*I,J)Y(J)=(1-TEMP)*FATHER(2*I-1,J)+TEMP*FATHER(2*I,J)WRITE(*,*)X(,J,)=,X(J)WRITE(*,*)Y(,J,)=,Y(J)CONTINUECALL GX1(YES(1) !调用约束函数(针对X)STOPWRITE(*,*)111CALL GX2(YES(2)!调用约束函数(针对Y)WRITE(*,*)YES(1)=,YES(1)WRITE(*,*)YES(2)=,YES(2)STOPWRITE(*,*)222IF (YES(1).EQ.1.AND.YES(2).EQ.1)THENDO 115 J=1 , NCHILD(2*I-1,J)= X(J)CHILD(2*I,J)=Y(J)CONTINUEELSEIF (YES(1).EQ.1)THENDO 120 J=1 , N7CHILD(2*I-1,J)= X(J)120125!110!CONTINUEELSEIF (YES(2).EQ.1)THENDO 125 J=1 , NCHILD(2*I,J)= Y(J)CONTINUEENDIFIF (YES(1).NE.1.OR.YES(2).NE.1)THENSTOPGOTO 130ENDIFCONTINUE !交叉操作结束STOP!注：交叉时父代应随机配对，如何实现？!现配对为1、2，3、4，（顺序配对）! DO 140 I=1 , 2*K! WRITE(5,*)I,(FATHER(I,J),J=1,N)!1050 FORMAT(1X,初始最优染色体为,3X,F10.6,3X,F10.6,3X,F10.6)!140 CONTINUE !输出交叉父代! WRITE(5,*)! DO 143 I=1 , 2*K! WRITE(5,*)I,(CHILD(I,J),J=1,N)!1060 FORMAT(1X,交叉操作CHILD(,I2,)=,3X,F10.6,3X,F10.6,3X,F10.6)!143 CONTINUE !输出交叉子代! WRITE(5,*)!以交叉子代代替V2中对应染色体DO 150 I=1,2*KDO 145 J=1,NV2(NUM(I),J)=CHILD(I,J)145150CONTINUECONTINUE!生成V3(实际上是将替换后的V2复制)DO 160 I =1,POPSIZEDO 155 J=1,NV3(I,J)=V2(I,J)155 CONTINUE! WRITE(5,1070)I,(V3(I,J),J=1,N)!输出V3!1070 FORMAT(1X,V3(,I2,)=,3X,F8.4,3X,F8.4,3X,F8.4)160!CONTINUE!计算各染色体的变异率DO 161 I=1,POPSIZEIF(FIT2(I).LT.FAVG2.OR.FMAX2.EQ.FAVG2)THENPM(I)=PM1ELSEPM(I)=PM1-(PM1-PM2)*(FIT2(I)-FAVG2)/(FMAX2-FAVG2)ENDIFWRITE(*,*)PM(,I,)=,PM(I)8161175180CONTINUE!计算各染色体的变异率结束!选择变异用染色体NUMBER=0DO 180 I=1,POPSIZETEMP=RAN3(IDUM)IF (TEMP.LT.PM(I) THENNUMBER=NUMBER+1NUM(NUMBER)=IDO 175 J=1 , NFATHER(NUMBER,J)=V2(I,J)CONTINUEENDIFCONTINUE !变异父代选择结束! DO 440 I=1 , NUMBER! WRITE(5,*)I,(FATHER(I,J),J=1,N)!5050 FORMAT(1X,变异操作FATHER(,I2,)=,3X,F10.6,3X,F10.6,3X,F10.6)!440 CONTINUE !输出变异父代! WRITE(5,*)!生成随机方向,变异时方向系数从0到M?185!190!CC210195DO 190 I=1, NTL(I)=L*1.0/GENWRITE(*,*)TL(I)CONTINUESTOP!变异操作#对基因上下限打折扣IF(NTIAO.EQ.1)THENRIGHT1(1)=(1-TL(1)/1.05)*RIGHTDIELSERIGHT1(1)=TL(1)*RIGHTDIIF(RIGHT1(1).LT.LEFT1(1)THENRIGHT1(1)=RIGHTDIENDIFENDIF#打折扣结束K=0DO 200 I=1,NUMBER!M=RAN3(IDUM)*MDO 195 J=1,NRLR=RAN3(IDUM)IF(RLR.GT.0.5)THEN !针对SIGN为0的情况YLR(J)=RIGHT1(J)-FATHER(I,J)X(J)=FATHER(I,J)+YLR(J)*(RLR*(1-TL(J)*BLRELSEYLR(J)=FATHER(I,J)-LEFT1(J) !针对SIGN为1的情况X(J)=FATHER(I,J)-YLR(J)*(RLR*(1-TL(J)*BLRENDIFCONTINUE9CALL GX1(YES(1) !调用约束函数!WRITE(*,*)YES(1)IF (YES(1).EQ.1)THEN!输出变异子代! WRITE(5,6060)I,(X(J),J=1,N)!6060 FORMAT(1X,变异操作CHILD(,I2,)=,3X,F8.4,3X,F8.4,3X,F8.4)!若变异后染色体满足约束条件，则直接替换V3中对应染色体DO 205 J=1 , NV3(NUM(I),J)=X(J)205!201200!248250255260CONTINUEELSEK=K+1IF (K.EQ.50)THEN!若同一方向上变异十次仍无满足约束条件染色体,则重新选择变异方向DO 201 J= 1 , ND(J)=RAN3(IDUM)CONTINUEK=0GOTO 210ELSEGOTO 210ENDIFENDIFCONTINUE!将新一代种群排序，计算其最优目标函数值,按从大到小排列DO 250 I=1 , POPSIZEWRITE(*,*)I=,INUM(I)=IDO 248 J=1 ,NX(J)=V3(I,J)WRITE(*,*)X(,J,)=,X(J)CONTINUECALL GX1(YES(I)!调用约束函数CALL MUBIAO(F(I) !调用目标函数CONTINUEDO 260 I=1 , POPSIZEMIN=F(I)NUMBER=NUM(I)DO 255 J=I , POPSIZEIF (F(J).LT.MIN) THENTEMP=MINMIN=F(J)F(J)=TEMPTEMP=NUMBERNUMBER=NUM(J)NUM(J)=TEMPENDIFCONTINUEF(I)=MINNUM(I)=NUMBERCONTINUEIF (BEST.GT.F(1) THEN!记录最小值10BEST=F(1)BEST1=LDO 271 I=1 , NBESTG(I)=V3(NUM(1),I)271275280!1080!290CONTINUEEND IFDO 280 I= 1,POPSIZE !(染色体排序后赋予V4形成新一代种群)DO 275 J=1 , NV4(I,J)=V3(NUM(I),J)CONTINUECONTINUEDO 290 I=1,POPSIZEWRITE(5,1080)I,(V4(I,J),J=1,N)FORMAT(1X,V4(,I2,)=,3X,F8.4,3X,F8.4,3X,F8.4)CONTINUE !输出V4!V4为新一代种群,将V4复制到VDO 300 I =1,POPSIZEDO 295 J=1,NV(I,J)=V4(I,J)295 CONTINUE300 CONTINUEDO I=1,NLMDBEST(I)=BESTG(I)! WRITE(*,*)LMDBEST(I)ENDDO!输出最优目标函数值! WRITE(5,1090) BEST!1090 FORMAT(1X,最优目标函数值=,3X,F8.6)WRITE(*,*)BEST=,BESTWRITE(*,*)BEST1=,BEST1WRITE(7,*)BEST1=,BEST1WRITE(5,*)(BESTG(J),J=1,N)!1100 FORMAT(1X,最优染色体为,3X,F10.6,3X,F10.6,3X,F10.6)! WRITE(8,*)遗传代数=,L! DO I=1,N! WRITE(8,*)I,LEFT1(I),RIGHT1(I)! ENDDO10000 CONTINUE !遗传循环结束! WRITE(*,*)自动画图!CALL HUAMIANTU! WRITE(*,*)122CLOSE(5)END !主程序结束C#生成滑动面的AUTOCAD脚本文件SUBROUTINE HUAMIANTUCOMMON/NB/NB !土坡表面几何特征点数NBCOMMON/NS/NS !土层数NSCOMMON/XU/XU(10000)!地表特征点水平坐标COMMON/YU/YU(1000,10000) !地表特征点下各土层界面的垂直坐标COMMON/WL/WL(10000) !地表特征点下地下水或浸润线的垂直坐标COMMON/XDING/XDING,YDING !读入坡顶的坐标值COMMON/N3/N3 !底部划分条块数11COMMON/N4/N4 !中部划分条块数COMMON/N5/N5 !顶部划分条块数COMMON/N6/N6 !第2种滑面顶部划分条块数1COMMON/N7/N7 !第2种滑面顶部划分条块数2COMMON/N1/N1 !土条划分数目COMMON/NTIAO/NTIAO !判断属于哪种分条模式COMMON/LMDBEST/LMDBEST(100)REAL XBEST(10000),YBEST(10000)REAL CIT3(10000),CIT5(10000)REAL CIT6(10000),CIT7(10000)REAL DX3(10000),DY3(10000)REAL DX6(10000),DY6(10000)REAL CIT3L,CIT5LREAL CIT6L,CIT7LREAL XJ3(10000),YJ3(10000)REAL XJ6(10000),YJ6(10000)REAL BL4,XL,YL !滑面下交点的坐标值REAL AA3,AA4 !地表特征线的斜率REAL LMDBESTPARAMETER(PII=3.1415926)IF(NTIAO.EQ.1)THENGOTO 361ELSEGOTO 461ENDIF361362363364365366N1=N3+N4+N5CIT3L=PII/(2*N3)CIT5L=PII/(2*N5)CIT3(1)=0.0DX3(1)=-1.0DY3(1)=0.0DO 362 I=2,N3+1CIT3(I)=CIT3(I-1)+CIT3LDX3(I)=-COS(CIT3(I)DY3(I)=-SIN(CIT3(I)CONTINUEDO 363 I=N3+2,N3+N4+1DX3(I)=0.0DY3(I)=-1.0CONTINUEDO 364 I=1,N3+N4+1CIT5(I)=0.0CONTINUEDO 365 I=N3+N4+2,N1+1CIT5(I)=CIT5(I-1)+CIT5LDX3(I)=SIN(CIT5(I)DY3(I)=-COS(CIT5(I)CONTINUEBL4=XDING/N4DO 366 I=1,N3+1XJ3(I)=0.0YJ3(I)=0.0CONTINUEJ=1DO 367 I=N3+2,N3+N4+112XJ3(I)=XJ3(I-1)+BL438863887!367!368!369!461!462463!IF(XJ3(I).LE.XU(J)GOTO 3887J=J+1IF(J.GT.NB)THENWRITE(*,*)找不到该点属于哪个特征段!STOPENDIFGOTO 3886AA3=(XJ3(I)-XU(J-1)/(XU(J)-XU(J-1)YJ3(I)=YU(1,J-1)+AA3*(YU(1,J)-YU(1,J-1)WRITE(*,*)YJ3(,I,)=,YJ3(I)CONTINUEDO 368 I=N3+N4+2,N1+1XJ3(I)=XDINGWRITE(*,*)XJ3(,I,)=,XJ3(I)YJ3(I)=YDINGWRITE(*,*)YJ3(,I,)=,YJ3(I)CONTINUERIGHT(1)=-XLXBEST(1)=-LMDBEST(1)YBEST(1)=0.0XBEST(N1+1)=LMDBEST(N1+1)+XDINGYBEST(N1+1)=YDINGXBEST(2)=XJ3(2)+DX3(2)*LMDBEST(2)YBEST(2)=YJ3(2)+DY3(2)*LMDBEST(2)DO 369 I=2,N1XBEST(I)=XJ3(I)+DX3(I)*LMDBEST(I)YBEST(I)=YJ3(I)+DY3(I)*LMDBEST(I)CONTINUEDO I=1,N1+1WRITE(6,*)XBEST(I),YBEST(I)ENDDOGOTO 4600ENDIF!过坡面的滑动面N1=N6+N7WRITE(*,*)N1=,N1WRITE(*,*)LMDBEST(1)J=1XL=XDING-LMDBEST(1)*XDING/SQRT(XDING*2+YDING*2)IF(XL.LE.XU(J)GOTO 463J=J+1IF(J.GT.NB)THENWRITE(*,*)找不到该点属于哪个特征段!STOPENDIFGOTO 462AA4=(XL-XU(J-1)/(XU(J)-XU(J-1)YL=YU(1,J-1)+AA4*(YU(1,J)-YU(1,J-1)WRITE(*,*)YJ3(,I,)=,YJ3(I)BET=ATAN(YDING-YL)/(XDING-XL)CIT6L=(PII/2.0-BET)/N6CIT7L=PII/(2.0*N7)13CIT6(1)=BETDO 464 I=2,N6+1CIT6(I)=CIT6(I-1)+CIT6LDX6(I)=-COS(CIT6(I)DY6(I)=-SIN(CIT6(I)464465466467!46846002900390032003600310033003400!3800!CONTINUEDO 465 I=1,N6+1CIT7(I)=0.0CONTINUEDO 466 I=N6+2,N1+1CIT7(I)=CIT7(I-1)+CIT7LDX6(I)=SIN(CIT7(I)DY6(I)=-COS(CIT7(I)CONTINUEDO 467 I=1,N1+1XJ6(I)=XDINGYJ6(I)=YDINGCONTINUEXBEST(1)=XLYBEST(1)=YLWRITE(*,*)XBEST(1),YBEST(1)XBEST(N1+1)=LMDBEST(N1+1)+XDINGYBEST(N1+1)=YDINGDO 468 I=2,N1XBEST(I)=XJ6(I)+DX6(I)*LMDBEST(I)YBEST(I)=YJ6(I)+DY6(I)*LMDBEST(I)CONTINUEDO I=1,N1+1WRITE(6,*)XBEST(I),YBEST(I)ENDDODO 3100 I=1,NSWRITE(3,2900)LINEFORMAT(5A)DO 3200 J=1,NBWRITE(3,3900)XU(J),YU(I,J)FORMAT(E11.5,1A,E11.5)CONTINUEWRITE(3,3600) FORMAT(1A)CONTINUEWRITE(3,2900)LINEDO 3300 I=1,NBWRITE(3,3900)XU(I),WL(I)CONTINUEWRITE(3,3600) WRITE(3,2900)LINEWRITE(3,3400)XBEST(1),YBEST(1)FORMAT(E11.5,1A,E11.5)WRITE(*,*