C BIOSA1 SOURCE CB215821 25/03/25 21:15:03 12217 SUBROUTINE BIOSA1(INDU,POTE,IMC,PN1,P1,P2,PP3,ROTA,VI,GY,IU,NB,TM, & UPSI,PENT1,PENT2,F1,F2,F3,ZF1,ZALF,ZDENS,ZBET,C, & ISECT,SANGLE,VX,VY,VZ,RESP) C Routine appelee par BIOSAV C Calcul du champ d'induction ou du potentiel vecteur en un point C cible PN1 IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C SANGLE est juste tansmis dans cette SUBROUTINE INTEGER SANGLE -INC PPARAM -INC CCOPTIO -INC CCREEL LOGICAL INDU,POTE DIMENSION C(3),PN1(3),RESP(3) DIMENSION VX(3),VZ(3),VY(3),VI(3),PP3(3) DIMENSION ZF1(1),ZALF(1),ZDENS(1),ZBET(1) DIMENSION P1(3),P2(3),P3(3),ROTA(3,3) DATA EPSIL3/1.D-12/,EPSIL2/1.D-8/ ALF=ZALF(1) BET=ZBET(1) DENS=ZDENS(1) XUN =1.D0 XDEUX =2.D0 XDEUX =2.D0 XDEMI =0.5D0 XPREC =1.D-6 C EMBRANCHEMENT VERS LA BONNE PARTIE DE LA SUBROUTINE C SELON LE TYPE D'INDUCTEUR (CERCLE, ARC, BARRE, FIL) GOTO (100,100,200,300),IMC C CAS DES CERCLES ET DES ARC 100 CONTINUE C MESURE DU POINT DANS LE REPERE LOCAL DO J=1,3 PN1(J)=PN1(J)-P1(J) ENDDO DO I= 1,3 VZ(I)=XZERO VY(I)=XZERO DO J= 1,3 VY(I)=VY(I)+ROTA(I,J)*PN1(J) VZ(I)=VZ(I)+ROTA(I,J)*PP3(J) ENDDO ENDDO C ON CHANGE DE RERERENTIEL C MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS ZP=VY(3) RP=SQRT(VY(1)**2+VY(2)**2) IF(A_EGALE_B (RP,XZERO)) THEN TMIN=XZERO COST=XUN SINT=XZERO RP=EPSIL2 ELSE TMIN=ACOS(VY(1)/(RP+EPSIL3)) COST=COS(TMIN) SINT=VY(2)/(RP+EPSIL3) ENDIF DO J= 1,3 C(J)=XZERO ENDDO C C CAS DU FIL CIRCULAIRE IF(BET.EQ.XZERO .AND. ALF.EQ.XUN) THEN RAYON=F1 COUR=DENS C ON CONSTRUIT LE REPERE LOCAL SIGT=XUN IF(ABS(SINT).GT.XPREC) SIGT=SIGN(XUN,-SINT) TMIN=SIGT*(TMIN-XPI)+XPI IF(IMC.EQ.1) TM=XDEUX*XPI TMAX=TMIN+TM CALL TEGRAL(SANGLE) IF(INDU) THEN CALL BFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF IF(POTE) THEN CALL AFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF RESP(1)=BX*COST-BY*SINT RESP(2)=BX*SINT+BY*COST RESP(3)=BZ C CAS DE LA PLAQUE CIRCULAIRE ELSE IF(BET.EQ.XZERO.AND.ISECT.EQ.3) THEN RINT=F1 REXT=F2 DR=REXT-RINT COUR=DENS*DR C ON CONSTRUIT LE REPERE LOCAL SIGT=XUN IF(ABS(SINT).GT.XPREC) SIGT=SIGN(XUN,-SINT) TMIN=SIGT*(TMIN-XPI)+XPI IF(IMC.EQ.1) TM=XDEUX*XPI TMAX=TMIN+TM CALL TEGRAL(SANGLE) IF(INDU) THEN CALL BPLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF IF(POTE) THEN CALL APLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF RESP(1)=BX*COST-BY*SINT RESP(2)=BX*SINT+BY*COST RESP(3)=BZ C CAS DU CYLINDRE ELSEIF(ALF.EQ.1D0) THEN RAYON=F1 H=F3 COUR=DENS*H C ON CONSTRUIT LE REPERE LOCAL SIGT=XUN IF(ABS(SINT).GT.XPREC) SIGT=SIGN(XUN,-SINT) TMIN=SIGT*(TMIN-XPI)+XPI IF(IMC.EQ.1) TM=XDEUX*XPI TMAX=TMIN+TM CALL TEGRAL(SANGLE) IF(INDU) THEN CALL BCYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF IF(POTE) THEN CALL ACYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF RESP(1)=BX*COST-BY*SINT RESP(2)=BX*SINT+BY*COST RESP(3)=BZ C AUTRES CAS ELSE IF (IMC.EQ.1.AND.ISECT.EQ.3) THEN IF(INDU) THEN CALL CHAMBO(IU,IC,NB,C(3),ZF1,ZALF,ZBET,ZDENS,ZP,RP,BZ,BR,BP) IF (RP.LE.UPSI) THEN RESP(1)=XZERO RESP(2)=XZERO ELSE RESP(1)=BR*VY(1)/RP RESP(2)=BR*VY(2)/RP ENDIF RESP(3)=BZ GOTO 999 ENDIF IF(POTE) THEN RINT=F1 REXT=F2 HCEN=F3 WIDTH=F2-F1 AIRE=F3*WIDTH COUR=DENS*AIRE DELTA=XZERO C ON CONSTRUIT LE REPERE LOCAL SIGT=XUN IF(ABS(SINT).GT.1.E-8) SIGT=SIGN(XUN,-SINT) TMIN=SIGT*(TMIN-XPI)+XPI TM=XDEUX*XPI TMAX=TMIN+TM PENT1=XZERO PENT2=XZERO CALL TEGRAL(SANGLE) CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT, & PENT1,PENT2,TMIN,TMAX,RP,ZP,BX,BY,BZ) RESP(1)=BX*COST-BY*SINT RESP(2)=BX*SINT+BY*COST RESP(3)=BZ GOTO 999 ENDIF ENDIF C CAS DES SECTIONS RECTANGULAIRES IF(ISECT.EQ.3) THEN C ON CALCULE POUR UN ARC ---> ANGLES A1=XZERO IF (ABS(VZ(1)).LT.UPSI) THEN A2=XPI/XDEUX ELSE AA=VZ(2)/VZ(1) A2=ATAN (AA) IF (A2.LT.XZERO) A2=XPI+A2 ENDIF IF(INDU) THEN CALL CHAMAR(IU,IC,C,A1,A2,F1,ALF,BET,DENS, & VY(1),VY(2),VY(3),BX,BY,BZ) RESP(1)=BX RESP(2)=BY RESP(3)=BZ ENDIF IF(POTE) THEN RINT=F1 REXT=F2 HCEN=F3 WIDTH=F2-F1 AIRE=F3*WIDTH COUR=DENS*AIRE PENT1=XZERO PENT2=XZERO C ON CONSTRUIT LE REPERE LOCAL SIGT=XUN IF(ABS(SINT).GT.XPREC) SIGT=SIGN(XUN,-SINT) TMIN=SIGT*(TMIN-XPI)+XPI TMAX=TMIN+TM CALL TEGRAL(SANGLE) CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT, & PENT1,PENT2,TMIN,TMAX,RP,ZP,BX,BY,BZ) RESP(1)=BX*COST-BY*SINT RESP(2)=BX*SINT+BY*COST RESP(3)=BZ ENDIF C CAS DES SECTIONS TRAPEZOIDALES ELSEIF(ISECT.EQ.1) THEN RINT=F1 REXT=F2 HCEN=F3 WIDTH=F2-F1 AIRE=F3*WIDTH COUR=DENS*AIRE DELTA=XDEMI*WIDTH*(PENT2-PENT1) HINT=HCEN-DELTA HEXT=HCEN+DELTA C ON CONSTRUIT LE REPERE LOCAL SIGT=XUN IF(ABS(SINT).GT.1.E-06) SIGT=SIGN(XUN,-SINT) TMIN=SIGT*(TMIN-XPI)+XPI IF(IMC.EQ.1) TM=XDEUX*XPI TMAX=TMIN+TM CALL TEGRAL(SANGLE) IF(HCEN.EQ.XZERO) THEN C CAS DU TRONC DE CONE TETA=ATAN(PENT1) COUR=DENS*WIDTH/COS(TETA) IF(INDU) THEN CALL BTRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX, & RP,ZP,BX,BY,BZ) ENDIF IF(POTE) THEN CALL ATRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX, & RP,ZP,BX,BY,BZ) ENDIF C AUTRES CAS ELSE IF(INDU) THEN CALL ARCTRA(SANGLE,COUR,HINT,HEXT,RINT,REXT, & PENT1,PENT2,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF IF(POTE) THEN CALL AARCTR(SANGLE,COUR,HINT,HEXT,RINT,REXT, & PENT1,PENT2,TMIN,TMAX,RP,ZP,BX,BY,BZ) ENDIF ENDIF RESP(1)=BX*COST-BY*SINT RESP(2)=BX*SINT+BY*COST RESP(3)=BZ ENDIF ENDIF IF(IERR.NE.0) RETURN GOTO 999 C CAS DES BARRES 200 CONTINUE C ON CHANGE DE RERERENTIEL DO I= 1,3 VY(I)=XZERO DO J= 1,3 VY(I)=VY(I)+ROTA(I,J)*PN1(J) ENDDO ENDDO C MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS C CAS DES SECTIONS RETANGULAIRES IF(ISECT.EQ.3) THEN IF(F2.EQ.XZERO) THEN C CAS DE LA PLAQUE RECTANGULAIRE DO I=1,3 VY(I)=XZERO DO J=1,3 VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J)) ENDDO ENDDO DY=F1 REXT=XDEMI*DY RINT=-REXT COUR=DENS*DY YMIN=-VY(1) YMAX=YMIN+GY XP=-VY(2) ZP=VY(3) IF(INDU) THEN CALL BPLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP) ENDIF IF(POTE) THEN CALL APLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP) ENDIF RESP(1)=BYP RESP(2)=-BXP RESP(3)=BZP ELSE IF(INDU) THEN CALL CHABAR(IU,VZ,GY,F1,F2,DENS,VY(1),VY(2),VY(3), & BXP,BYP,BZP) RESP(1)=BXP RESP(2)=BYP RESP(3)=BZP ENDIF IF(POTE) THEN DO I=1,3 VY(I)=XZERO DO J=1,3 VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J)) ENDDO ENDDO DY=F1 DZ=F2 REXT=XDEMI*DY RINT=-REXT AIRE=DY*DZ COUR=DENS*AIRE PENT1=XZERO PENT2=XZERO YMIN=-VY(1) YMAX=YMIN+GY XP=-VY(2) ZP=VY(3) CALL ASEGTR(COUR,DZ,DZ,RINT,REXT,PENT1,PENT2, & YMIN,YMAX,XP,ZP,BXP,BYP,BZP) RESP(1)=BYP RESP(2)=-BXP RESP(3)=BZP ENDIF ENDIF C CAS DES SECTIONS TRAPEZOIDALES ELSEIF(ISECT.EQ.1) THEN DO I=1,3 VY(I)=XZERO DO J=1,3 VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J)) ENDDO ENDDO DY=F1 DZ=F2 REXT=XDEMI*DY RINT=-REXT AIRE=DY*DZ COUR=DENS*AIRE DELTA=XDEMI*DY*(PENT2-PENT1) HINT=DZ-DELTA HEXT=DZ+DELTA YMIN=-VY(1) YMAX=YMIN+GY XP=-VY(2) ZP=VY(3) IF(INDU) THEN CALL SEGTRA(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2, & YMIN,YMAX,XP,ZP,BXP,BYP,BZP) ENDIF IF(POTE) THEN CALL ASEGTR(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2, & YMIN,YMAX,XP,ZP,BXP,BYP,BZP) ENDIF RESP(1)=BYP RESP(2)=-BXP RESP(3)=BZP ENDIF GOTO 999 C CAS DU FIL 300 CONTINUE GY=XZERO DO I=1,3 VY(I)=P2(I)-P1(I) VI(I)=PN1(I)-P1(I) GY=GY+(VY(I)*VY(I)) ENDDO GY=SQRT(GY) CALL PROVEC(VI,VY,VZ) CALL PROVEC(VY,VZ,VX) CALL ROTREP (ROTA,VX,VY,VZ) IF(IERR.NE.0) RETURN XP=XZERO YMIN=XZERO DO J= 1,3 XP=XP+ROTA(1,J)*VI(J) YMIN=YMIN-ROTA(2,J)*VI(J) ENDDO YMAX=YMIN+GY COUR=DENS IF(INDU) THEN CALL BFILDR(COUR,YMIN,YMAX,XP,BZP) RESP(1)=XZERO RESP(2)=XZERO RESP(3)=BZP ENDIF IF(POTE) THEN CALL AFILDR(COUR,YMIN,YMAX,XP,AYP) RESP(1)=XZERO RESP(2)=AYP RESP(3)=XZERO ENDIF 999 CONTINUE RETURN END