(*
 * Do Not Remove This Notice
* Per Unit C and L for Microstrip Lines.
* Converted from Fortran to Mathematica by
* Sasan Ardalan, December 11, 2018
* http://www.radiocalc.com  *)


(* From Fortran Code MSTRP.FOR by Clayton Paul:

To Determine the NxN Capacitance and Inductance Matrices
  of a N-Conductor Coupled Microstrip Consisting of N Strips of
  Equal-Width above an Infinite Ground Plane.The Number
  of Divisions per Strip for the Charge Distribution
 on each Land is Denoted by NCDIV.The Maximum Problem Size is
  MSIZE=N*NCDIV.C Output Data are in PUL.DAT and Input Data are in MSTRP.IN
 
  Written by:Clayton R.Paul
  Department of Electrical Engineering
  University of Kentucky
 Lexington,KY 40506
 
*)



Phi[D_, DELTA_, FAC_,ER_, T_]:=
Module[{xx=0,DP,DM,TOW,PHI},
   DP=2.*D/DELTA+1.;
       DM=2.*D/DELTA-1.;
       TOW=4.*T/DELTA;
If[D== 0,
PHI=FAC*(0.5*Log[ 1.+TOW^2]+TOW*ArcTan[1./TOW]),

  PHI=DM*Log[DM]-DP*Log[DP]-0.5*DM*Log[DM^2+TOW^2]+0.5*DP*Log[DP^2+TOW^2]+TOW*(ArcTan[DP/TOW]-ArcTan[DM/TOW]);
        PHI=PHI*FAC/2 
        
 
]


]  (* End Phi    *)


PhiBB[D_, DELTA_, FAC_,ER_, T_]:=
Module[{xx=0,DP,DM,TOW,XK,EREFF,F1,F2,TNT,NN,BEFORE,TEMP,P0,PhiVal},
   DP=2.*D/DELTA+1.;
       DM=2.*D/DELTA-1.;
       TOW=4.*T/DELTA;
XK=(ER-1.)/(ER+1.);
EREFF=(ER+1)/2.;

 

If[D== 0,
BEFORE=0.;
NN=0;
NN=NN+1;
TNT=TOW*NN;
F1=XK^(NN-1);
F1=F1*((-1.)^NN);
F2=0.5*Log[1.+TNT^2]+(TNT)*ArcTan[1./TNT];
TEMP=F1*F2;
           BEFORE=BEFORE+TEMP;
 

(* 200 *)

While[(Abs[TEMP/BEFORE]  >= 10^(-5)), 
NN=NN+1;
TNT=TOW*NN;
F1=XK^(NN-1);
F1=F1*((-1.)^NN);
F2=0.5*Log[1.+TNT^2]+(TNT)*ArcTan[1./TNT];
TEMP=F1*F2;
               BEFORE=BEFORE+TEMP   
];

(*
IF(ABS(TEMP/BEFORE).GE.1.E-5) GO TO 200
*)
PhiVal=-BEFORE*FAC*(1.+XK)/EREFF
,
(*  ELSE *)


    P0=0.5*DM*Log[DM]-0.5*DP*Log[DP];
            BEFORE=0.;
             NN=0;
TEMP=0;
NN=NN+1;
TNT=TOW*NN;
F1= XK^(NN-1);
F1=F1*((-1.0)^NN);
F2=0.25*DM*Log[DM^2+TNT^2]-0.25*DP*Log[DP^2+TNT^2]+0.5*TNT*(ArcTan[DM/TNT]-ArcTan[DP/TNT]);
                 TEMP=F1*F2;
                 BEFORE=BEFORE+TEMP;
 
While[(Abs[TEMP/BEFORE]  ≥ 10^(-5)) ,
(*100 Continue *)
NN=NN+1;
TNT=TOW*NN;
F1= XK^(NN-1);
F1=F1*((-1.0)^NN);
F2=0.25*DM*Log[DM^2+TNT^2]-0.25*DP*Log[DP^2+TNT^2]+0.5*TNT*(ArcTan[DM/TNT]-ArcTan[DP/TNT]);
                 TEMP=F1*F2;
                 BEFORE=BEFORE+TEMP
] ;
(* IF(ABS(TEMP/BEFORE).GE.1.E-5) GO TO 100  *)

PhiVal=P0+BEFORE*(1.+XK);
PhiVal=PhiVal*FAC/EREFF
]
]   (* End PhiBB   *)


(*********************************************************)
(* START PROGRAM                                          *)


NCDIV=50;  (*NUMBER OF CONDUCTOR SUBSECTIONS *)
(* Setup Input      ****************)
(*********************************************************)
(* Specify Muticonductor Microstrip    System             *)

NN=2 ; (* TOTAL NUMBER OF LANDS (EXCLUSIVE OF GND PLANE) *)
 W=189; (* LAND WIDTH (mils) *)
SEP=189.; (* EDGE-TO-EDGE LAND SEPARATION (mils) *)
T=61.;        (*   BOARD THICKNESS (mils)  *)
ER=2.2 ;       (* BOARD RELATIVE DIELECTRIC CONSTANT *)
(*************** End Specification       ******************)
(*********************************************************)


(*********** DO NOT EDIT BELOW THIS LINE                     *****)
CMTM=2.54*10^(-5);

W=W*CMTM;
SEP=SEP*CMTM;

T=T*CMTM;


(******* END  Setup Input      ****************)
V0=2.997925*10^8;
V02=V0*V0;
EPS=1.0/(V02*4*10^(-7)*Pi);
  
FAC=N[1./(2.*Pi*EPS)];





DELTA=W/NCDIV;


NTOT=NN*NCDIV;
(*    Test 
Z1=Phi[0, DELTA, FAC,ER, T]
Z2=PhiBB[0, DELTA, FAC,ER, T]
Print["Z1=",Z1," Z2=", Z2]
*)
Clear[AMat]
Clear[CapMat]
Clear[InductMat]
Clear[Cap0Mat]




MSIZE=150; (* Parameter *)
AMat=Array[ #1+#2&,{NTOT,NTOT}];
CapMat=Array[ #1+#2&,{NN,NN}];
InductMat=Array[ #1+#2&,{NTOT,NTOT}];
Cap0Mat=Array[ #1+#2&,{NN,NN}];



(******************WITHOUT DIELECTRIC*****************************************FILL A1**********************************************************)
DIST=0;
For[ii=1,ii<NCDIV+1,ii++,
          AMat[[1,ii]]=Phi[DIST, DELTA, FAC,ER, T];
 
         AMat[[ii,1]]=  AMat[[1,ii]];
         DIST=DIST+DELTA
]
For[ii=2,ii<NCDIV+1,ii++,
For[jj=ii,jj<NCDIV+1,jj++,
AMat[[ii,jj]]=AMat[[1,jj-ii+1]];
AMat[[jj,ii]]=AMat[[ii,jj]]
]
]
 

(******FILL A2--AN******************************************************)
For [ii=2,ii<NN+1, ii++,
DIST=(ii-1)*(W+SEP);
For[jj=1, jj<NCDIV+1,jj++,
XJMI=jj-1;
AMat[[1,jj+(ii-1)*NCDIV]]=Phi[DIST+XJMI*DELTA, DELTA, FAC,ER, T];
 ];
For[jj=2, jj<NCDIV+1,jj++,
XJMI=jj-1;
AMat[[jj, (ii-1)*NCDIV+1]]=Phi[DIST-XJMI*DELTA, DELTA, FAC,ER, T] 
 
];
For[kk=2, kk<NCDIV+1,kk++,
For[mm=kk, mm<NCDIV+1,mm++,
AMat[[kk,mm+(ii-1)*NCDIV]]=AMat[[(kk-1),(mm-1)+(ii-1)*NCDIV]];
AMat[[mm,kk+(ii-1)*NCDIV]]=AMat[[(mm-1),(kk-1)+(ii-1)*NCDIV]]
]
]
]
 
(*******FILL AT2-ATN*****************************************************)

For[ii=1, ii<NCDIV+1,ii++,
For[jj=ii, jj<NTOT+1,jj++,
AMat[[jj,ii]]=AMat[[ii,jj]]
];
]
 
(********FILL A MATRIX****************************************************)
NM1=NN-1;
For[ii=1, ii<NM1+1,ii++,
NMI=NN-ii;
For[jj=1, jj<NMI+1,jj++,
For[kk=1, kk<NCDIV+1,kk++,
For[ll=1, ll<NCDIV+1,ll++,
AMat[[kk+jj*NCDIV,ll+jj*NCDIV+(ii-1)*NCDIV]]=AMat[[kk,ll+(ii-1)*NCDIV]]
];
];
];  
]
NM2=NN-2;
 
For[ii=1, ii<NM2+1,ii++,
NMI=NN-(ii+1);

Print [NMI];
For[jj=1, jj<NMI+1,jj++,
For[kk=1, kk<NCDIV+1,kk++,
For[ll=1, ll<NCDIV+1,ll++,
AMat[[kk+(jj+ii)*NCDIV,ll+jj*NCDIV]]=AMat[[kk+ii*NCDIV,ll]]
          ];
             ];
                ]  


]
 
AMat= Inverse[AMat];
 
 

For[ii=1, ii<NN+1,ii++,
For[jj=1, jj<NN+1,jj++,
SUM=0;
For[kk=1, kk<NCDIV+1,kk++,
For[mm=1, mm<NCDIV+1,mm++,
   SUM=SUM+AMat[[kk+(ii-1)*NCDIV,mm+(jj-1)*NCDIV]]


 
];
];
 
Cap0Mat[[ii,jj]]=SUM;


];
]
 
CapMat=Cap0Mat
 
CapMat=Inverse[CapMat];
V02Inv=1/V02;
 
InductanceMatrix=V02Inv*CapMat;
Print["Inductance Matrix"]
Print[InductanceMatrix]


 
(******************WITH DIELECTRIC**************************************
******FILL A1**********************************************************)
DIST=0;
For[ii=1, ii<NCDIV+1,ii++,
   AMat[[1,ii]]=PhiBB[DIST, DELTA, FAC,ER, T];
   AMat[[ii,1]]=AMat[[1,ii]];

     
DIST=DIST+DELTA
] 
 
For[ii=2, ii<NCDIV+1,ii++,
For[jj=ii, jj<NCDIV+1,jj++,
AMat[[ii,jj]]=AMat[[1,jj-ii+1]];
AMat[[jj,ii]]=AMat[[ii,jj]]
];
]
  
(*******FILL A2--AN******************************************************)
For[ii=2, ii<NN+1,ii++,
DIST=(ii-1)*(W+SEP);
                      For[jj=1, jj<NCDIV+1,jj++,
                             XJM1=jj-1;
                            AMat[[1,jj+(ii-1)*NCDIV]]=PhiBB[DIST+XJM1*DELTA, DELTA, FAC,ER, T]
                  ];
                     For[jj=2, jj<NCDIV+1,jj++,
                             XJM1=jj-1;
                            AMat[[jj,(ii-1)*NCDIV+1]]=PhiBB[DIST-XJM1*DELTA, DELTA, FAC,ER, T]
                  ];
                     For[kk=2, kk<NCDIV+1,kk++,
                         For[mm=kk, mm<NCDIV+1,mm++,
                           AMat[[kk,mm+(ii-1)*NCDIV]]=AMat[[(kk-1),(mm-1)+(ii-1)*NCDIV]];
                           AMat[[mm,kk+(ii-1)*NCDIV]]=AMat[[mm-1,(kk-1)+(ii-1)*NCDIV]]
                     ];
                  ];
]
  
(*******FILL AT2-ATN*****************************************************)
NTOT=NN*NCDIV;
For[ii=1, ii<NCDIV+1,ii++,
     For[jj=1, jj<NTOT+1,jj++,
       AMat[[jj,ii]]=AMat[[ii,jj]]
];
]
   
(********FILL A MATRIX****************************************************)
NM1=NN-1;
For[ii=1, ii<NM1+1,ii++,
NMI=NN-ii;
For[jj=1, jj<NMI+1,jj++,
For[kk=1, kk<NCDIV+1,kk++,
For[ll=1, ll<NCDIV+1,ll++,
AMat[[kk+jj*NCDIV,ll+jj*NCDIV+(ii-1)*NCDIV]]=AMat[[kk,ll+(ii-1)*NCDIV]]
];
];
];
]
 
NM2=NN-2;
For[ii=1, ii<NM2+1,ii++,
NMI=NN-(ii+1);
For[jj=1, jj<NMI+1,jj++,
For[kk=1, kk<NCDIV+1,kk++,
For[ll=1, ll<NCDIV+1,ll++,
AMat[[kk+(jj+ii)*NCDIV,ll+jj*NCDIV]]=AMat[[kk+ii*NCDIV,ll]]
];
];
];
]

  

(*******DETERMINE CAPACITANCE MATRIX WITH DIELECTRIC**********************)
AMat=Inverse[AMat];
 

For[ii=1, ii<NN+1,ii++,
For[jj=1, jj<NN+1,jj++,
SUM=0;
For[kk=1, kk<NCDIV+1,kk++,
For[mm=1, mm<NCDIV+1,mm++,
   SUM=SUM+AMat[[kk+(ii-1)*NCDIV,mm+(jj-1)*NCDIV]]
];
];
 
CapMat[[ii,jj]]=SUM;

];
]
Print["Capacitance Matrix"]
Print[CapMat]
Print["Capacitance (Cap0)  Matrix"] 
Print[Cap0Mat]