eigen value and eigen vector

bekar__bahadır1 · ‎02-05-2019

I'm working on quantum structures, so I need to calculate eigenvalues and eigenvalues in a very precise drawing, I want to work with values such as n = 10000 and above. Mkl library that I use it is not very favorable to offer an alternative program. 4 core 3.6 ghz 12gb ram on laptop,
I also have a 200-core 120 gb ram li host computer.

my code

c use imsl
INCLUDE 'LINK_FNL_STATIC.H'
USE EVESB_INT
C USE EVESF_INT
c USE EPISF_INT
c USE EPISB_INT
C USE CSDER_INT
C USE CSINT_INT

    IMPLICIT NONE
INTEGER I,II,J,K,L,M,N,NDATA,NINTV,LDA,LDEVEC,NCODA,NEVAL,
$NEVEC,INT_TIME,IK,MXEVAL
PARAMETER (M=101,N=M*M,NCODA=M,NDATA=N,LDA=N,LDEVEC=N
$,NEVEC=4,MXEVAL=4)
C=======================================================================
   REAL*8 A(LDA,N),ALPHA,AALPHA,PI,BREAK(NDATA),AA,BB,OTOP1U,SEBIN,
$XX,YY,ZZ,EPSILON,LAMDA,DLAMDA,RRO,VB(N),PSO(N),VDC(M,M),C,VO,
$VVO,DU,RO1RO2,RO3,RO4,A1,A2,A3,A4,B1,B2,B3,B4,PSU(M,M),
$PSIN(M,M),PSF(M,M),DRO,TOP,RO,EM,VL(N),DZ,AYIL,RYIL,F,ETA,B,GAMA,
$RR,U,H,XI,YI,KZZ,INTEN,HPLANCK,VVVO,DX,DY,X(M),Y(M),MY,RRIC,RIC,KZ
$,P,EO,EB,EIK,EUC,VM(N),TOPKISI,KISIBIR,KISIUC,EPS,FXSU,FXOU,OTOPXU
$,TOPXU,BETA1,BETA3,TOPBETA,TOPYU,M12,HW,T,EF,EIN,BETA3U,BETA3A,BET
$A1U,BETA1A,E,FXOA1,FXSA1,FXOA2,FXSA2,NR,R,OTOPXA1,OTOPXA2,TOPXA1,T
$OPXA2,TOPYA1,TOPYA2,TOPSON,EPSO,EVAL(NEVEC),EVEC(LDEVEC,NEVEC),F1
$U,TOP1U,F1SU,OSI,VS(N),EYUKU,SIGMA,TZAMAN,CISIK,KISIBIRA,KISIBIR
$U,KISIUCU,KISIUCA,RR1,R1,RR2,R2,TTB,TB,PII,XLAMDA,DXLAMDA,SAY,OTOP
$XA,TOPYA,OEBIN,FK0,FK1,TOPXA,FXSA,ATA,FXOA,PS,RO1,RO2,FXSA3,
$TOPXA3,OTOPXA3,OTOP1OU,OTOP2OU,F2U,F2SU,OTOP2U,M11,M22,KS1U,
$KS31A,KS31U,KS32U,KS33U,KS34U,KS32A,KS33A,TOPKS,TOP2U,BETA31U,
$BETA32U,BETA33U,BETA31A,KS1A,KS1,KS3
   REAL*8 OTOPYU,OTOPYA,TOPZU,TOPZA,Z,LL,LA,OTOPYA3,TOPYA3,FXOA3
C=======================================================================
LOGICAL SMALL
CHARACTER*8 CHAR_TIME
CALL TIME(CHAR_TIME)
WRITE(*,*)'TIME1=', CHAR_TIME
   PI=4.0D0*DATAN(1.0D0)
C=======================================================================
C OPEN(1,FILE='ALGAAS k L60.DAT')
c   OPEN(2,FILE='ALGAAS s-PISI_L20.DAT')
c OPEN(3,FILE='ALGAAS silindir A-E.DAT')
c   OPEN(4,FILE='ALGAAS PISILER DELTOID00t70.DAT')
c   OPEN(5,FILE='AlGaAs r2 B VE ENERJI M1 KZ2.DAT')
C   OPEN(6,FILE='GAALAS IC BARIYER VE TABAN ENERJI.DAT')
c   OPEN(7,FILE='GAALAS BAGLANMA ENERJI-TBB LAZER 00KARE XIYI00.DAT')
   OPEN(8,FILE='AlGaAs 1-2 BETA INTEN03 ALFA 60 k.DAT')
   OPEN(9,FILE='AlGaAs 1-2 K_INDEX INTEN03 ALFA 60 k.DAT')

C%%%%%%%%%%%%%%%%%%%%%%%%%%%% SABITLER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
INTEN=0.300000E10 !MEGA WATT/ CMKARE----> METREYEKAREYE CEVRILIYORc
AALPHA=60.0D0 !LAZER GENLIĞI (ANG)
LL=105.0
C=================================ALGAAS===================================
MY=0.06650D0
   EPSILON=10.9 ! 13.18 10.90D0
C****************** ALGAN************************************************
C MY=0.13D0
C EPSILON=9.7 !YÜKSEK FREKANS 5.3!STATIK9.7
C===================GaINAS ===============================================
c MY=0.023+0.037*0.3+0.003*(0.3)**2
C EPSILON=15.1-2.87*0.3+0.67*(0.3)**2!STATIK
c EPSILON=12.3-1.4*0.3 !YUKSEK
C==========================================================================
   RYIL=(13605.698110D0*MY/(EPSILON**2))
   AYIL=0.52917724820D0*EPSILON/MY

C==================POT=================================================
c VO=228000000000000.0D0/RYIL !ALGAAS IçIN
VO=228.0/RYIL !DİKKAT
C VO=345.0D0/RYIL !ALGAN IçIN
c VO=227.70D0/RYIL !GAINAS IçIN
c DO 5555 TTB=5.0,200,5.0
c   DO 5555 AALPHA=0,100,5.0
c DO 5555 B=0.0D0,20.0D0,1.0D0
C DO 5555 RRIC=50.0,150.0,5.0
C DO 5555 VVVO=100.0,300.0,10.0

F=00.0D0 !ELEKTRIK ALAN ŞIDDETI (KV/CM)
   B=0.0D0 !MANYETIK ALAN ŞIDDETI (TESLA)

   XI=0.0000001D0 !YAB. AT. KONUMU
   YI=0.0000001D0 !YAB. AT. KONUMU
   RR=220.0012345670D0 !Dış GENIşLIK
   RRIC=50.00D0 !Iç KUYU GENIşLIğI
TTB=50.0
C======================= OPTIK GEçIş KATSAYıLARı ======================
EYUKU=1.60217733E-19!DSQRT(2.0D0) !COULOMB
SIGMA=3.0E22 !M-^3 TAŞIYICI YOĞUNLUĞU için

TZAMAN=5E12 !PIKO SANIYE SANIYEYE CEVRILIP CARPIM DURUMUNDA
C TZAMAN=(1.0/1.5)*1E12 !Algan
CISIK=2.99792458E8 !METRE/SANIYE
EPSO=8.854187817E-12!C^2/(NEWTON.METREKARE)
NR=3.2!DSQRT(EPSILON)
HPLANCK=1.05457266E-34!J.SANIYE
C****************** DONUSUMLER **************************************
C**********************************************************************
   ALPHA=AALPHA/AYIL
LA=LL/AYIL
   RIC=RRIC/AYIL
   R=RR/AYIL
RR1=40.0
RR2=150.0
R1=RR1/AYIL
R2=RR2/AYIL
TB=TTB/AYIL
C VO=VVVO/RYIL
KZ=0.0D0/(Ric)!DALGA SAYıSı
   EM=0.0D0 !AZIMUTHAL MAGNETIK ALAN
C======================================================================
   ETA=0.010D0*AYIL*F/RYIL
GAMA=4.254381195E-6*EPSILON*EPSILON*B/(MY*MY)
   DX=(2.0D0*R)/REAL(M-1)
   DY=(2.0D0*R)/REAL(M-1)
    DZ=(2.0D0*R)/REAL(M-1)
DRO=R/REAL(M-1)
   AA=4.0D0/(DX*DX)
   BB=-1.0D0/(DX*DX)
C***********************************************************************
   PRINT*,'EPSILON:',EPSILON,'EPSILON=10,89 ISE LAZER AKTIF'
PRINT*,'RYIL:',RYIL
PRINT*,'MY:',MY
   PRINT*,'AALPHA',AALPHA,'ANGUSTRON'
PRINT*,'B=',B
PRINT*,'M=',EM
PRINT*,'KZ=',KZ
PRINT*,'RRIC',RRIC
PRINT*,'GAMMA=',GAMA
PRINT*,'VO=',VO*RYIL
PRINT*,'INTEN=',INTEN

C************** AZIMUTHAL MAGNETIK ALAN *******************************
   II=1
DO K=1,M
   X(K)=-R+REAL(K-1)*DX
IF(ABS(X(K)).LE.0.000000001) GO TO 32
32 DO L=1,M
Y(L)=-R+REAL(L-1)*DY
IF(ABS(Y(L)).LE.0.0000000001) GOTO 33
RO=DSQRT(X(K)*X(K)+Y(L)*Y(L))
IF(RO.LT.RIC)THEN
VB(II) =((EM*EM/(RO*RO))+KZ*KZ-(GAMA*RO*RO*KZ/RIC)
$+(GAMA*GAMA*(RO**4)/(4*(RIC**2))))
VS(II)=0.0!0.25*GAMA*GAMA*RO*RO
ELSE
VB(II)=0.0!(EM*EM/(RO*RO))+KZ*KZ+2*KZ*GAMA*RIC*LOG(RIC/RO)+
C $GAMA*GAMA*RO*RO*(LOG(RIC/RO)**2)
VS(II)=0.0!0.25*GAMA*GAMA*RO*RO
END IF
33 II=II+1
   END DO
END DO
C********************* LAZER GIYDIRILIYOR ******************************
   PRINT*,'LAZER GIYDIRILIYOR'
   II=1
   DO 999 L=1,M
YY=Y(L)
DO 888 K=1,M
   XX=X(K)
C=======================================================================
   TOP=0.0D0
   DU=0.0010D0
   DO U=0.0D0,2.0D0*PI,DU
TOPSON=VVO(XX+ALPHA*DSIN(U),YY,RYIL,AYIL,R1,R2,TB,RIC,LA,VO)
   TOP=TOP+(TOPSON)*DU
   END DO
   TOP=TOP/(2.0D0*PI)
C=======================================================
   VDC(K,L)=TOP
    VL(II)=VDC(K,L)
VM(II)=VL(II)+VB(II)+VS(II)
II=II+1
WRITE(1,19)X(K)*AYIL,Y(L)*AYIL,VDC(K,L)*RYIL
888   CONTINUE
999   CONTINUE
    PRINT*, 'KUYU TANIMLANDI', 'VB'

19 FORMAT(3(2X,F14.8))
18 FORMAT(5(2X,F14.8))
177 FORMAT(4(2X,F14.8))
C********************* MATRIS *****************************************
A=0.0D0 !AMATRIS BOLOK(2M+1,N=M*M)
DO L=M+1,M*M
A(1,L)=BB
ENDDO

DO L=2,M*M
   A(M,L)=BB
ENDDO
   DO L=M+1,M*M-1,M
   A(M,L)=0.0D0
ENDDO

               DO L=1,M*M
               A(M+1,L)=AA+Vm(L)
               ENDDO
C                       DO L=1,M*M-1
C                       A(M+2,L)=BB
C                       END DO
C                               DO L=1,M*M-M
C                               A(2*M+1,L)=BB
C                               ENDDO
C*********** MATRIS EKRANA YAZDıRıLıYOR *******************************
c DO 40 K=1,M+1
c WRITE(*,'(1X,6(F6.1))') (A(K,L),L=1,m)
c WRITE(1,'(1X,6(F6.1))') (A(K,L),L=1,m)
c40 CONTINUE
c PAUSE
c   STOP
C***********************************************************************
SMALL =.TRUE.
CALL DEVESB(N,NEVEC,A,LDA,NCODA,SMALL,EVAL,EVEC,LDEVEC)
C CALL DEVESF (N, NEVEC, A, LDA, SMALL, EVAL, EVEC, LDEVEC)
C PRINT*,'PERFORMANS INDEX=',PII
C PII=EPISB(NEVEC,A,NCODA,EVAL,EVEC)
c PII= EPISF(NEVEC,A,EVAL,EVEC)
CALL TIME(CHAR_TIME)
WRITE(*,*)'TIME2=', CHAR_TIME

   EO=(EVAL(NEVEC))*RYIL !TABAN DURUM MEV CINSINDEN
   EB=(EVAL(NEVEC-1))*RYIL !1.UYARıLMıS DURUM MEV CINSINDEN
   EIK=(EVAL(NEVEC-2))*RYIL !2.UYARıLMıS DURUM MEV CINSINDEN
EUC=(EVAL(NEVEC-3))*RYIL !3.UYARıLMıS DURUM MEV CINSINDEN
C=======================================================================
   II=1
   DO L=1,M
   DO K=1,M
   PSIN(K,L)=EVEC(II,NEVEC)*AYIL*1E-10 !PISILER METRE BOYUTUNDA
   PSF(K,L)=EVEC(II,NEVEC-1)*AYIL*1E-10 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PSU(K,L)=EVEC(II,NEVEC-2)*AYIL*1E-10
PSO(II)=EVEC(II,NEVEC)
WRITE(4,511)X(K)*AYIL,Y(L)*AYIL,PSIN(K,L)/1E-10,PSF(K,L)/1E-10,
$EVEC(II,NEVEC-2)*AYIL*1E-10/1E-10
   II=II+1
ENDDO
   ENDDO
   PRINT*, 'EO=',EO,'MEV'
   PRINT*, 'E1=',EB,'MEV'
   PRINT*, 'E2=',EIK,'MEV'
PRINT*, 'E3=',EUC,'MEV'
   WRITE(3,188)AALPHA,EO,EB,EIK,EUC
c WRITE(5,188)B,EO,EB,EIK,EUC
c WRITE(5,*)B,EO
C WRITE(6,*)RRIC,EO
   EF=EB
   EIN=EO

C********************* M12 HESABI ********************************
OTOPXU=0.0D0
OTOP1OU=0.0D0
OTOP2OU=0.0D0
   OTOPXA1=0.0D0
   OTOPXA2=0.0D0
OTOPXA3=0.0
   TOPYU=0.0D0
TOP1U=0.0D0
TOP2U=0.0D0
   TOPYA1=0.0D0
   TOPYA2=0.0D0
TOPYA3=0.0
   DO 400 L=1,M
   YY=Y(L)
C===================INTEGRALIN X KıSMı BAşLıYOR======================
   FXOU=0.0D0
F1U=0.0D0
F2U=0.0D0

   FXOA1=0.0D0
   FXOA2=0.0D0
FXOA3=0.0
       TOPXU=0.0D0
TOP1U=0.0D0
TOP2U=0.0D0
   TOPXA1=0.0D0
TOPXA2=0.0D0
TOPXA3=0.0D0
II=1
RO=DSQRT(XX*XX+YY*YY)

DO 300 K=1,M
   XX=X(K)*AYIL*1E-10 !METREYE CEVIRDIK
yy=y(K)*AYIL*1E-10
FXSU=PSIN(K,L)*xx*PSF(K,L) !X YONDE POLRIZASYON ıKEN XX Y ıSE YY KULLAN
F1SU=PSIN(K,L)*xx*PSIN(K,L)
F2SU=PSF(K,L)*xx*PSF(K,L)

   FXSA1=PSIN(K,L)*PSIN(K,L)
   FXSA2=PSF(K,L)*PSF(K,L)
FXSA3=PSU(K,L)*PSU(K,L)

   TOPXU=TOPXU+(FXSU+FXOU)*(DX)/2.0D0
TOP1U=TOP1U+(F1SU+F1U)*DX/2.0D0
TOP2U=TOP2U+(F2SU+F2U)*DX/2.0D0
   TOPXA1=TOPXA1+(FXSA1+FXOA1)*(DX)/2.0D0
   TOPXA2=TOPXA2+(FXSA2+FXOA2)*(DX)/2.0D0
TOPXA3=TOPXA3+(FXSA3+FXOA3)*(DX)/2.0D0
   FXOU=FXSU
F1U=F1SU
F2U=F2SU
   FXOA1=FXSA1
   FXOA2=FXSA2
FXOA3=FXSA3
300 CONTINUE
C==================== X KıSMı BITTI ==============================
TOPYU=TOPYU+(OTOPXU+TOPXU)*(DY)/2.0D0
TOP1U=TOP1U+(OTOP1U+TOP1U)*DY/2.0D0
TOP2U=TOP2U+(OTOP2U+TOP2U)*DY/2.0D0
       TOPYA1=TOPYA1+(OTOPXA1+TOPXA1)*(DY)/2.0D0
   TOPYA2=TOPYA2+(OTOPXA2+TOPXA2)*(DY)/2.0D0
TOPYA3=TOPYA3+(OTOPXA3+TOPXA3)*(DY)/2.0D0
   OTOPXU=TOPXU
OTOP1U=TOP1U
OTOP2U=TOP2U
       OTOPXA1=TOPXA1
   OTOPXA2=TOPXA2
OTOPXA3=TOPXA3
400 CONTINUE
C%%%%%%%%%%%%%%% NORMALIZE PISI %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   DO K=1,M
   DO L=1,M
   WRITE(2,17)X(K)*AYIL,Y(L)*AYIL,(PSIN(K,L)/DSQRT(TOPYA1))**2,
$(PSF(K,L)/DSQRT(TOPYA2))**2,(PSU(K,L)/DSQRT(TOPYA3))**2,VDC(K,L)
$*RYIL
   ENDDO
   ENDDO
PRINT*, 'TOPYA3',TOPYA3
PRINT*, 'TOPYA2',TOPYA2
PRINT*, 'TOPYA1',TOPYA1
M12=(TOPYU*EYUKU)/(DSQRT(TOPYA1)*DSQRT(TOPYA2))
M11=(TOPYU*EYUKU)/(DSQRT(TOPYA1)*DSQRT(TOPYA1))
M22=(TOPYU*EYUKU)/(DSQRT(TOPYA2)*DSQRT(TOPYA2))

c OSI=(2*MY*9.1093897E-31/(HPLANCK**2))*((EB-EO)*1.6021773E-22)*
c $(((TOPROU)/(DSQRT((TOPYA1))*DSQRT((TOPYA2))))**2)
c   PRINT*, 'M12=',M12
c WRITE(13,*)B,OSI
C WRITE(10,*)B,VS(II)
C PRINT *, 'B=',B, OSI
C PAUSE
C STOP
C////////////////////// OPTIK GECIS ///////////////////////////////
C########### 1. VE 3. DERECE ABSORTSIYON KATSAYISI ####################

C######################################################################
C########### 1. VE 3. DERECE ABSORTSIYON KATSAYISI ####################
DO 4444 HW=0.0D0,250.0D0,1.0D0
BETA1U=SIGMA*(HW*1.6021773E-22)*(M12*M12)*TZAMAN
BETA1A=CISIK*EPSO*NR*(((((EF-EIN-HW)*1.6021773E-22)**2))+
$(HPLANCK*TZAMAN)**2)

BETA1=BETA1U/BETA1A

BETA3U=INTEN*2.0*SIGMA*((M12)**4)*(HW*1.6021773E-22)*TZAMAN
BETA31U=DABS((M22-M11)/(2.0*M12))**2
BETA32U=(((EF-EIN-HW)*1.6021773E-22)**2)-(HPLANCK*TZAMAN)**2
BETA33U=2*((EF-EIN)*1.6021773E-22)*((EF-EIN-HW)*1.6021773E-22)

BETA3A=CISIK*CISIK*EPSO*EPSO*NR*NR*((((((EF-EIN-HW)*1.6021773E-22)
$)**2)+(HPLANCK*TZAMAN)**2)**2)
BETA31A=((EF-EIN)*1.6021773E-22)**2+(HPLANCK*TZAMAN)**2


BETA3=-(BETA3U/BETA3A)*(1-(((BETA31U)*(BETA32U+BETA33U))/BETA31A))

   TOPBETA= BETA1+BETA3

C***************************************************************
C*************** DIREK ABSORTSION KATSAYISI ********************
   KS1U=SIGMA*((EF-EIN-HW)*1.6021773E-22)*((M12)**2)
KS1A=2.0*NR*NR*EPSO*((((EF-EIN-HW)*1.6021773E-22)**2
$)+(HPLANCK*TZAMAN)**2)

   KS1=KS1U/KS1A

KS31U=INTEN*SIGMA*((EF-EIN-HW)*1.6021773E-22)*((M12)**4)
KS31A=NR*NR*NR*EPSO*EPSO*CISIK*((((EF-EIN-HW)*1.6021773E-22)**2
$+(HPLANCK*TZAMAN)**2)**2)

KS32U=DABS((M22-M11)/(2.0*M12))**2
KS33U=((EF-EIN)*1.6021773E-22)*((EF-EIN-HW)*1.6021773E-22)**2
KS34U=((HPLANCK*TZAMAN)**2)*(3*((EF-EIN)*1.6021773E-22)-2*
$(HW*1.6021773E-22))
KS32A=((EF-EIN)*1.6021773E-22)**2+(HPLANCK*TZAMAN)**2
KS33A=(EF-EIN-HW)*1.6021773E-22

KS3=-(KS31U/KS31A)*(1-KS32U*((KS33U-KS34U)/(KS32A*KS33A)))

   TOPKS=KS1+KS3

c   PRINT*,'TOPKISI',TOPKS
WRITE(8,51)HW,BETA1/1e4,BETA3/1e4,TOPBETA/1e4
   WRITE(9,51)HW,Ks1,Ks3,TOPKS
c WRITE(12,*)INTEN/1E7,TOPBETA/1E2
4444 CONTINUE !HW FOTON ENERJISI DöNGUSU   , ALPHA DONGUSU
C################## YABANCI ATOM #########################################

C XLAMDA=0.1
C   DXLAMDA=0.010D0
C   SAY=0.
C   OEBIN=-1.0D30
C150 CONTINUE
C=====================================================integralin Z kısmı başlıyor==========
C OTOPYU=0.0D0
C   OTOPYA=0.0D0
C
C   TOPZU=0.0D0
C   TOPZA=0.0D0
C DZ=0.10D0

C       DO 500 Z=-R,R,DZ
C IF(ABS(Z).LE.0.000000010D0)GOTO 500
C
C
C
C
C=====================================================integralin y kısmı başlıyor==========
C OTOPXU=0.0D0
C   OTOPXA=0.0D0
C   TOPYU=0.0D0
C TOPYA=0.0D0
C II=1
C
C DO 400 J=1,M
C YY=Y(J)
C====================================================integralin x kısmı başlıyor===========
C   FXOU=0.
C   FXOA=0.
C   TOPXU=0.
C   TOPXA=0.
C
C
C   DO 300 I=1,M
C   XX=X(I)
C   RO1=DSQRT((XX-XI+ALPHA)**2+(YY-YI)**2+Z*Z)
C RO2=DSQRT((XX-XI-ALPHA)**2+(YY-YI)**2+Z*Z)
C           PS=EVEC(II,NEVEC)*DEXP(-DABS(RO1)+DABS(RO2)/(2.0D0)*XLAMDA)

C ATA=((1.0D0/RO1)+(1.0D0/RO2))/2.0D0
C FXSU=(PS*ATA*PS)
C   FXSA=(PS*PS)

C   TOPXU=TOPXU+(FXSU+FXOU)*DX/2.0D0
C TOPXA=TOPXA+(FXSA+FXOA)*DX/2.0D0

C   II=II+1

C   FXOU=FXSU
C   FXOA=FXSA
C WRITE(*,*)'RO1',RO1
C WRITE(*,*)'RO2',RO2
C WRITE(*,*)'PS',PS
C WRITE(*,*)'ata',ata
C WRITE(*,301)Z,XX,YY,TOPXA,TOPXU
C300 CONTINUE
C301 FORMAT(5(2X,F10.6))
C=================================================== x kısmı bitti=========================

C TOPYU=TOPYU+(OTOPXU+TOPXU)*DY/2.0D0
C   TOPYA=TOPYA+(OTOPXA+TOPXA)*DY/2.0D0

C   OTOPXU=TOPXU
C   OTOPXA=TOPXA

C400 CONTINUE
C=================================================== y kısmı bitti=========================

C TOPZU=TOPZU+(OTOPYU+TOPYU)*DZ/2.0D0
C   TOPZA=TOPZA+(OTOPYA+TOPYA)*DZ/2.0D0

C   OTOPYU=TOPYU
C   OTOPYA=TOPYA

C500 CONTINUE
C=================================================== Z kısmı bitti=========================

C SEBIN=-(1.0D0/XLAMDA**2.)+2.0D0*(TOPZU/TOPZA) !bağlanma enerjisiC
C WRITE(*,*)XLAMDA,SEBIN,SAY

C PAUSE
C STOP

C====================================bağlanma enerjisi için hassaslaştırma yapılıyor=======
C IF(SEBIN.LT.OEBIN)THEN
C       IF(SAY.GT.5)GO TO 250
C           DXLAMDA=-DXLAMDA/5.0D0
C           SAY=SAY+1
C   ENDIF

C XLAMDA=XLAMDA+DXLAMDA
C   OEBIN=SEBIN

C   GO TO 150

C250 CONTINUE
C===========================================bağlanma enerjisi daha hassas bulundu=========
C========================================
C CALL TIME(char_time)
C WRITE(*,*)'TIME3=', char_time
C WRITE(7,*)ttb,SEBIN*RYIL
C WRITE(*,*)TTB,SEBIN*RYIL

C========================================

C700 CONTINUE

51 FORMAT(4(1X,F15.11))
511 FORMAT(5(2X,F25.19))
17 FORMAT(6(2X,F20.14))
16 FORMAT(3(2X,F14.8))
188 FORMAT(5(2X,F14.8))
C PAUSE
C STOP
c5555 CONTINUE
PAUSE
   STOP
   END
C=======================================================
C============================ FUNCTIONS ================
C=======================================================
C=======================================================
   FUNCTION VVO(XX,YY,RYIL,AYIL,R1,R2,TB,RIC,LA,VO)
   IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 LA
c==================== deltoid bariyerli ==================
c rdis=150.0/ayil
c if (abs(xx).ge.(abs(rdis)-abs(yy)))vvo=vo
c if(abs(xx).lt.(abs(rdis)-abs(yy)).and.abs(xx).gt.(abs(ric+tb)-
c $abs(yy)))vvo=0.0
c if(abs(xx).le.(abs(ric+tb)-abs(yy)).and.abs(xx).ge.(abs(ric)-
c $abs(yy)))vvo=vo
c if(abs(xx).lt.(abs(ric)-abs(yy)))VVO=0.0

C============== KARE =====================================
IF (ABS(XX).Le.(LA/2.0).AND.ABS(YY).Le.(LA/2.0))THEN
VVO=0.0
ELSE
VVO=VO
END IF

C******************** UCGEN*****************************************
c IF(ABS(XX).LT.(100/AYIL))THEN
c IF(ABS(yy).LE.ABS((50./AYIL)+(xx/2.0)))THEN
c VVO=0.0
c ELSE
c VVO=VO
c END IF
c ELSE
c VVO=VO
c END IF

c////////////////////DELTOİD/////////////////////////////////////////
c IF(ABS(YY).LT.(LA/SQRT(2.0)))THEN
c IF(ABS(XX).GT.ABS((LA/SQRT(2.0))-ABS(YY)))THEN
c VVO=VO
c ELSE
c VVO=0.0
c END IF
c ELSE
c VVO=VO
c END IF

C333333333333333333333333333333 DOUBLE KARE 333333333333333333333333333333333333333
C IF(ABS(YY).LT.(75./AYIL))THEN
C IF(ABS(XX).LT.(150.0/AYIL).AND.ABS(XX).GT.(50./AYIL))THEN
C VVO=0.0
C ELSE
C VVO=VO
C END IF
C ELSE
C VVO=VO
C ENDIF



C22222222222222222222222222222222 SILINDIR barıyerli2222222222222222222222222
C   RO=SQRT((XX)**2+(YY)**2)
C IF(RO.GE.R2) VVO=VO
C IF(RO.LT.R2.AND.RO.GT.(R1+TB))VVO=0.0
C IF(RO.LE.(R1+TB).AND.RO.GE.R1)VVO=VO
C IF(RO.LT.R1) VVO=0.0
C111111111111111111111 1SİLİNDİR 1111111111111111111111111111111
c    RO=SQRT((XX)**2+(YY)**2)
c    IF(RO.GT.(LA/2.0))THEN
c VVO=VO
c ELSE
c VVO=0.0D0
c ENDIF
C000000000000000000000 PARABOLIC 00000000000000000000000000000
C   IF(RRO.GE.RIC)THEN
C               VVO=VO
C   ELSE
C               VVO=VO*((1/RIC)**2)*RRO*RRO
C               ENDIF
C0000000000000000000 PETEK 00000000000000000000000000000000000000000000
C VVO=VO
C A1=150.0/AYIL !X1 KOORDINATI
C A2=-150.0/AYIL!X2 KOORDINATI
C A3=-150.0/AYIL!X3 KOORDINATI
C A4=150.0/AYIL!X4 KOORDINATI
C B1=150.0/AYIL!Y1 KOORDINATI
C B2=150.0/AYIL!Y2 KOORDINATI
C B3=-150.0/AYIL!Y3 KOORDINATI
C B4=-150.0/AYIL!Y4 KOORDINATI
C R1=25.0/AYIL!1 NOLU DAIRE Y CAPI
C R2=25.0/AYIL!2 NOLU DAIRE Y CAPI
C R3=25.0/AYIL!3 NOLU DAIRE Y CAPI
C R4=25.0/AYIL!4 NOLU DAIRE Y CAPI
C R5=50.0/AYIL!4 NOLU DAIRE Y CAPI
C
C RO=SQRT((XX)**2+(YY)**2)
C RO1=SQRT((XX-A1)**2+(YY-B1)**2)
C RO2=SQRT((XX-A2)**2+(YY-B2)**2)
C RO3=SQRT((XX-A3)**2+(YY-B3)**2)
C RO4=SQRT((XX-A4)**2+(YY-B4)**2)
C IF(RO1.LT.R1)VVO=0.0
C IF(RO2.LT.R2)VVO=0.0
C IF(RO3.LT.R3)VVO=0.0
C IF(RO4.LT.R4)VVO=0.0
C IF(RO.LT.r5)VVO=0.0


RETURN
   END

bekar__bahadır1 · ‎02-05-2019