      SUBROUTINE TYPE701(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: TYPE701

C  COPYRIGHT 2007 NATIONAL RENEWABLE ENERGY LABORATORY
C  Greg Barker's Model of the Sandia formula
C  Photovoltaic Module using Sandia model.


C INPUT #


C 1  IcBeam           (kJ/hr m2)
C 2  IcDiffuse        (kJ/hr m2)
C 3  IncidenceAngle   (deg)
C 4  ZenithAngle      (deg)
C 5  Tamb             (deg C)
C 6  WindSpeed        (m/s)
C 7  Altitude         (m)
C 8  OperatingVoltage (Volts)


C PARAMETER #

C 1  NumSeries        (unitless) (number of modules in series in one string)
C 2  NumStrings       (unitless) (number of strings in parallel)
C 3  Acoll            (m2, single module)
C 4  NcellSer         (unitless)
C 5  NparSerCells     (unitless)
C 6  Isc0             (Amps)
C 7  Voc0             (Volts)
C 8  Imp0             (Amps)
C 9  Vmp0             (Volts)
C 10 aIsc             (1/degC)
C 11 aImp             (1/degC)
C 12 c_0              (unitless)
C 13 c_1              (unitless)
C 14 BVoc0            (Volts/degC)
C 15 mBVoc            (Volts/degC)
C 16 BVmp0            (Volts/degC)
C 17 mBVmp            (Volts/degC)
C 18 DiodeFactor      (unitless)
C 19 c_2              (unitless)
C 20 c_3              (unitless)
C 21 a_0              (unitless)
C 22 a_1              (unitless)
C 23 a_2              (unitless)
C 24 a_3              (unitless)
C 25 a_4              (unitless)
C 26 b_0              (unitless)
C 27 b_1              (unitless)
C 28 b_2              (unitless)
C 29 b_3              (unitless)
C 30 b_4              (unitless)
C 31 b_5              (unitless)
C 32 DT0              (deg C)
C 33 fd               (unitless)
C 34 a                (unitless)
C 35 b                (unitless)
C 36 c_4              (unitless)
C 37 c_5              (unitless)
C 38 Ix0              (Amps)
C 39 Ixx0             (Amps)
C 40 c_6              (unitless)
C 41 c_7              (unitless)
C 42 IrradianceThreshold (kJ/hr m2) [number above which a fit to the single-diode model will be attempted]
C 43 InverterCapacityFactor (unitless)


C OUTPUT #

C 1  OperatingVoltage (Volts)
C 2  OperatingCurrent (Amps)
C 3  OperatingPower   (kJ/hr)
C 4  Eff              (unitless)
C 5  Vmp              (Volts)
C 6  Imp              (Amps)
C 7  Pmp              (kJ/hr)
C 8  EffMax           (unitless)
C 9  Isc              (Amps)
C 10 Voc              (Volts)
C 11 Tcell            (deg C)
C 12 Tback            (deg C)
C 13 AMa              (unitless)
C 14 F1               (unitless)
C 15 F2               (unitless)
C 16 Ix               (Amps)
C 17 Vx               (Volts)
C 18 Ixx              (Amps)
C 19 Vxx              (Volts)
C 20 Converged        (1=yes,0=no)
C 21 MaxError         (fraction)
C 22 Rs_init          (ohms)
C 23 z_init           (unitless)
C 24 TotIter          (unitless)
C 25 Rs               (ohms)
C 26 z                (unitless)
C 27 Rsh              (ohms)
C 28 Io               (amps)
C 29 IL               (amps)
C 30 AbsErrorVx       (kJ/hr)
C 31 AbsErrorVmp      (kJ/hr)
C 32 AbsErrorVxx      (kJ/hr)
C 33 FracErrorVx      (kJ/hr)
C 34 FracErrorVmp     (kJ/hr)
C 35 FracErrorVxx     (kJ/hr)
C 36 Minimum Inverter Capacity(kJ/hr)

C  ACCESS TRNSYS FUNCTIONS
	USE TrnsysFunctions
	USE TrnsysConstants

C STANDARD TRNSYS DECLARATIONS      
      IMPLICIT NONE						!REQUIRES THE USER TO DEFINE ALL VARIABLES BEFORE USING THEM
	DOUBLE PRECISION XIN				!THE ARRAY FROM WHICH THE INPUTS TO THIS TYPE WILL BE RETRIEVED
	DOUBLE PRECISION OUT				!THE ARRAY WHICH WILL BE USED TO STORE THE OUTPUTS FROM THIS TYPE
	DOUBLE PRECISION TIME				!THE CURRENT SIMULATION TIME - YOU MAY USE THIS VARIABLE BUT DO NOT SET IT!
	DOUBLE PRECISION PAR				!THE ARRAY FROM WHICH THE PARAMETERS FOR THIS TYPE WILL BE RETRIEVED
	DOUBLE PRECISION T					!AN ARRAY CONTAINING THE RESULTS FROM THE DIFFERENTIAL EQUATION SOLVER
	DOUBLE PRECISION DTDT				!AN ARRAY CONTAINING THE DERIVATIVES TO BE PASSED TO THE DIFF.EQ. SOLVER
	DOUBLE PRECISION TIME0,TFINAL,DELT	!SIMULATION START TIME, STOP TIME, AND TIMESTEP
	INTEGER*4 ICNTRL					!AN ARRAY FOR HOLDING VALUES OF CONTROL FUNCTIONS WITH THE NEW SOLVER
      INTEGER*4 INFO(15)					!THE INFO ARRAY STORES AND PASSES VALUABLE INFORMATION TO AND FROM THIS TYPE
	INTEGER*4 NUMPAR,NUMIN, NUMOUT      !Number of Parameters, Inputs, and Outputs
	INTEGER*4 NoutToCheck               !Number of Outputs for checking connections using
	INTEGER*4 NoutActual                !Number of Outputs in the OUT() array
      CHARACTER*3 OCHECK					!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE OUTPUTS
	CHARACTER*3 YCHECK					!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE INPUTS

      PARAMETER (NUMIN=8,NUMPAR=43,NUMOUT=36)

      DIMENSION XIN(NUMIN)
      DIMENSION OUT(NUMOUT)
	DIMENSION PAR(NUMPAR)
	DIMENSION YCHECK(NUMIN),OCHECK(NUMOUT)

C Intermediate calculation variables
	INTEGER IUNIT
      INTEGER MaxIter
      DOUBLE PRECISION Ee
      DOUBLE PRECISION FitTolerance
      DOUBLE PRECISION Rs_Step
      DOUBLE PRECISION z_Step
      DOUBLE PRECISION Rs_InitialMatrixWidth
      DOUBLE PRECISION z_InitialMatrixWidth
      DOUBLE PRECISION TRW_Ix
      DOUBLE PRECISION TRW_Ixx
	DOUBLE PRECISION SingleModuleVoltage
	DOUBLE PRECISION MinIrradiance
	DOUBLE PRECISION I_TRW
	DOUBLE PRECISION Error
	DOUBLE PRECISION AbsError
	DOUBLE PRECISION Current_x
	DOUBLE PRECISION Current_mp
	DOUBLE PRECISION Current_xx
	DOUBLE PRECISION MaxTRWerror
	DOUBLE PRECISION RDCONV
	DOUBLE PRECISION PI
	DOUBLE PRECISION LUW
	LOGICAL UseSingleDiodeModel

C INPUTS      
	DOUBLE PRECISION IcBeam                   !(kJ/hr m2)
      DOUBLE PRECISION IcDiffuse, IcDiffuseCPV  !(kJ/hr m2)
      DOUBLE PRECISION IncidenceAngle           !(deg)
      DOUBLE PRECISION ZenithAngle      !(deg)
      DOUBLE PRECISION Tamb             !(deg C)
      DOUBLE PRECISION WindSpeed        !(m/s)
      DOUBLE PRECISION Altitude         !(m)
      DOUBLE PRECISION OperatingVoltage !(Volts)
      
C PARAMETERS
      DOUBLE PRECISION NumSeries        !(unitless)
	DOUBLE PRECISION NumStrings       !(unitless)
      DOUBLE PRECISION Acoll            !(m2)
      DOUBLE PRECISION NcellSer         !(unitless)
      DOUBLE PRECISION NparSerCells     !(unitless)
      DOUBLE PRECISION Isc0             !(Amps)
      DOUBLE PRECISION Voc0             !(Volts)
      DOUBLE PRECISION Imp0             !(Amps)
      DOUBLE PRECISION Vmp0             !(Volts)
      DOUBLE PRECISION aIsc             !(1/degC)
      DOUBLE PRECISION aImp             !(1/degC)
      DOUBLE PRECISION c_0              !(unitless)
      DOUBLE PRECISION c_1              !(unitless)
      DOUBLE PRECISION BVoc0            !(Volts/degC)
      DOUBLE PRECISION mBVoc            !(Volts/degC)
      DOUBLE PRECISION BVmp0            !(Volts/degC)
      DOUBLE PRECISION mBVmp            !(Volts/degC)
      DOUBLE PRECISION DiodeFactor      !(unitless)
      DOUBLE PRECISION c_2              !(unitless)
      DOUBLE PRECISION c_3              !(unitless)
      DOUBLE PRECISION a_0              !(unitless)
      DOUBLE PRECISION a_1              !(unitless)
      DOUBLE PRECISION a_2              !(unitless)
      DOUBLE PRECISION a_3              !(unitless)
      DOUBLE PRECISION a_4              !(unitless)
      DOUBLE PRECISION b_0              !(unitless)
      DOUBLE PRECISION b_1              !(unitless)
      DOUBLE PRECISION b_2              !(unitless)
      DOUBLE PRECISION b_3              !(unitless)
      DOUBLE PRECISION b_4              !(unitless)
      DOUBLE PRECISION b_5              !(unitless)
      DOUBLE PRECISION DT0              !(deg C)
      DOUBLE PRECISION fd               !(unitless)
      DOUBLE PRECISION a                !(unitless)
      DOUBLE PRECISION b                !(unitless)
      DOUBLE PRECISION c_4              !(unitless)
      DOUBLE PRECISION c_5              !(unitless)
      DOUBLE PRECISION Ix0              !(Amps)
      DOUBLE PRECISION Ixx0             !(Amps)
      DOUBLE PRECISION c_6              !(unitless)
      DOUBLE PRECISION c_7              !(unitless)
	DOUBLE PRECISION IrradianceThreshold !(kJ/hr m2)
	DOUBLE PRECISION InverterCapacityFactor !(unitless)

C OUTPUTS
C	DOUBLE PRECISION OperatingVoltage !(Volts) (Already Declared)
	DOUBLE PRECISION OperatingCurrent !(Amps)
	DOUBLE PRECISION OperatingPower   !(kJ/hr)
      DOUBLE PRECISION Eff              !(unitless)
	DOUBLE PRECISION Vmp              !(Volts)
	DOUBLE PRECISION Imp              !(Amps)
	DOUBLE PRECISION Pmp              !(kJ/hr)
      DOUBLE PRECISION EffMax           !(unitless)
	DOUBLE PRECISION Isc              !(Amps)
	DOUBLE PRECISION Voc              !(Volts)
	DOUBLE PRECISION Tcell            !(deg C)
	DOUBLE PRECISION Tback            !(deg C)
	DOUBLE PRECISION AMa              !(unitless)
	DOUBLE PRECISION F1               !(unitless)
	DOUBLE PRECISION F2               !(unitless)
	DOUBLE PRECISION Ix               !(Amps)
	DOUBLE PRECISION Vx               !(Volts)
	DOUBLE PRECISION Ixx              !(Amps)
	DOUBLE PRECISION Vxx              !(Volts)
      LOGICAL Converged     !(1=yes,0=no)
      DOUBLE PRECISION MaxError         !(fraction)
      DOUBLE PRECISION Rs_init          !(ohms)
      DOUBLE PRECISION z_init           !(unitless)
      INTEGER TotIter       !(unitless)
      DOUBLE PRECISION Rs               !(ohms)
      DOUBLE PRECISION z                !(unitless)
      DOUBLE PRECISION Rsh              !(ohms)
      DOUBLE PRECISION Io               !(amps)
      DOUBLE PRECISION IL               !(amps)
      DOUBLE PRECISION AbsErrorVx       !(kJ/hr)
      DOUBLE PRECISION AbsErrorVmp      !(kJ/hr)
      DOUBLE PRECISION AbsErrorVxx      !(kJ/hr)
      DOUBLE PRECISION FracErrorVx      !(kJ/hr)
      DOUBLE PRECISION FracErrorVmp     !(kJ/hr)
      DOUBLE PRECISION FracErrorVxx     !(kJ/hr)
	DOUBLE PRECISION MinInverterCapacity !(kJ/hr)
      DOUBLE PRECISION dTcorr, Tref
C FUNCTIONS
      DOUBLE PRECISION TRW_Current
	DOUBLE PRECISION SandiaModuleTemperature
	DOUBLE PRECISION SandiaTcellFromTmodule
	DOUBLE PRECISION AbsoluteAirMass
	DOUBLE PRECISION SandiaF1
	DOUBLE PRECISION SandiaF2
	DOUBLE PRECISION SandiaIsc
	DOUBLE PRECISION SandiaEffectiveIrradiance
	DOUBLE PRECISION SandiaImp
	DOUBLE PRECISION SandiaVoc
	DOUBLE PRECISION SandiaVmp
	DOUBLE PRECISION SandiaIx
	DOUBLE PRECISION SandiaIxx
	DOUBLE PRECISION SingleDiodeCurrentTRW
      DOUBLE PRECISION EffInterpolate
	
CAD   storage for multiple point efficiency changes (2009.9.14)	
      DOUBLE PRECISION MpeRadArr(5),MpeEffArr(5), TotalRad

      DATA RDCONV/0.017453/,PI/3.1415927/

C SET THE VARIABLE TYPES FOR THE INPUTS AND OUTPUTS
C (NOTE: VOLTS AND AMPS ARE SET AS "DIMENSIONLESS=DM1"

C UNITS FOR INPUTS:
      DATA YCHECK/
     1 'IR1','IR1','DG1','DG1','TE1','VE1','LE1','DM1'/

C UNITS FOR OUTPUTS:
      DATA OCHECK/
	1 'DM1','DM1','PW1','DM1','DM1','DM1','PW1','DM1','DM1','DM1',
     1 'TE1','TE1','DM1','DM1','DM1','DM1','DM1','DM1','DM1','DM1',
     1 'DM1','DM1','DM1','DM1','DM1','DM1','DM1','DM1','DM1','PW1',
	1 'PW1','PW1','DM1','DM1','DM1','PW1'
	1/

C-----------------------------------------------------------------------------------------------------------------------
C    SET THE VERSION INFORMATION FOR TRNSYS
      IF(INFO(7).EQ.-2) THEN
	   INFO(12)=16
	   RETURN 1
	ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C FIRST CALL OF SIMULATION - DO ALL THE BOOKKEEPING WORK (ONLY ONCE)
      IF(INFO(7).EQ.-1) THEN

C SET THE UNIT NUMBER SO IN FUTURE ITERATIONS, THE READING OF THE 
C PARAMETERS MAY BE SKIPPED TO SAVE TIME
        IUNIT=0

C INFO(6) CONTAINS THE NUMBER OF OUTPUTS
	  INFO(6)=NUMOUT

C CALL THE TYPECK ROUTINE TO MAKE SURE THE INPUT FILE HAS THE CORRECT NUMBER
C OF INPUTS, PARAMETERS, AND DERIVATIVES.
        CALL TYPECK(1,INFO,NUMIN,NUMPAR,0)

C CALL THE RCHECK ROUTINE TO CHECK UNIT TYPES OF INPUTS AND OUTPUTS
        CALL RCHECK(INFO,YCHECK,OCHECK)

      ENDIF !FIRST CALL OF SIMULATION

C SEE IF THE PARAMETERS NEED TO BE READ AT THIS ITERATION AND
C SET THE CURRENT UNIT NUMBER IF THEY DO
      IF(INFO(1) .NE. IUNIT) THEN
  
        IUNIT=INFO(1)

C GET THE PARAMETERS
        NumSeries        =DBLE(PAR(1))
	  NumStrings       =DBLE(PAR(2))
	  Acoll            =DBLE(PAR(3))
        NcellSer         =DBLE(PAR(4))
        NparSerCells     =DBLE(PAR(5))
        Isc0             =DBLE(PAR(6))
        Voc0             =DBLE(PAR(7))
        Imp0             =DBLE(PAR(8))
        Vmp0             =DBLE(PAR(9))
        aIsc             =DBLE(PAR(10))
        aImp             =DBLE(PAR(11))
        c_0              =DBLE(PAR(12))
        c_1              =DBLE(PAR(13))
        BVoc0            =DBLE(PAR(14))
        mBVoc            =DBLE(PAR(15))
        BVmp0            =DBLE(PAR(16))
        mBVmp            =DBLE(PAR(17))
        DiodeFactor      =DBLE(PAR(18))
        c_2              =DBLE(PAR(19))
        c_3              =DBLE(PAR(20))
        a_0              =DBLE(PAR(21))
	  a_1              =DBLE(PAR(22))
        a_2              =DBLE(PAR(23))
        a_3              =DBLE(PAR(24))
        a_4              =DBLE(PAR(25))
        b_0              =DBLE(PAR(26))
        b_1              =DBLE(PAR(27))
        b_2              =DBLE(PAR(28))
        b_3              =DBLE(PAR(29))
        b_4              =DBLE(PAR(30))
        b_5              =DBLE(PAR(31))
        DT0              =DBLE(PAR(32))
        fd               =DBLE(PAR(33))
        a                =DBLE(PAR(34))
        b                =DBLE(PAR(35))
        c_4              =DBLE(PAR(36))
        c_5              =DBLE(PAR(37))
        Ix0              =DBLE(PAR(38))
        Ixx0             =DBLE(PAR(39))
        c_6              =DBLE(PAR(40))
        c_7              =DBLE(PAR(41))
	  IrradianceThreshold = DBLE(PAR(42))
	  InverterCapacityFactor=DBLE(PAR(43))
      
	ENDIF !PARAMETERS NEED TO BE READ IN

C GET THE INPUTS AT THE CURRENT ITERATION AND TIMESTEP
      IcBeam           =XIN(1) !(kJ/hr m2)
      IcDiffuse        =XIN(2) !(kJ/hr m2)
      IncidenceAngle   =XIN(3) !(deg)
      ZenithAngle      =XIN(4) !(deg)
      Tamb             =XIN(5) !(deg C)
      WindSpeed        =XIN(6) !(m/s)
      Altitude         =XIN(7) !(m)
      OperatingVoltage =XIN(8) !(Volts)

      MinIrradiance    =1.0 !Assume no operation if Ic is below this number (kJ/hr m2)
      
      
C  MpeRadArr is in Watt/m2
        MpeRadArr(1) = b_0
        MpeRadArr(2) = b_1
        MpeRadArr(3) = b_2
        MpeRadArr(4) = b_3
        MpeRadArr(5) = b_4

C  MpeEffArr is in fraction (0..1)        
        MpeEffArr(1) = c_0
        MpeEffArr(2) = c_1
        MpeEffArr(3) = c_2
        MpeEffArr(4) = c_3
        MpeEffArr(5) = c_4

CNB 4-16-06 The following code (IF statement) is for flat-plate 1pt-eff PV (CPV below)
CNB 12-23-04  Single-efficiency model option put in to allow simple analysis with SAM
      IF ((a_0.LT.0).AND.(a_3.GE.0.)) Then
CNB    a_0 is the negative fractional efficiency of the PV panel
CC Calculate max efficiency:
         EffMax = -a_0
                  
C   Convert radiation into W/m2 for interpolation (1 Watt = 3.6 kJ/hr)        
         TotalRad = (IcBeam+IcDiffuse)/3.6
         
C   Linearly interpolate efficiency from radiation level         
         EffMax = EffInterpolate( TotalRad, MpeRadArr, MpeEffArr )
                  
C  Pmp is the only real output available from this simple model
CNB This is only for a single module         Pmp = EffMax*(IcBeam+IcDiffuse)*Acoll
         Pmp = (EffMax*(IcBeam+IcDiffuse)*AColl)*NumStrings*NumSeries

        If (Pmp > 0.001) THEN
CNB 7-25-05 Adding code to correct for temperature correction
CNB use King method to calculate cell temperature
CNB Note that need appropriate fd, a and b for this call
C Calculate back-of-module temperature:
          Tback = SandiaModuleTemperature(IcBeam/3.6,IcDiffuse/3.6,
	1                                WindSpeed,Tamb,fd,a,b)

CNB Note that need fd and DT0 for this call
C Calculate cell temperature:
          Tcell = SandiaTcellFromTmodule(Tback,IcBeam/3.6,IcDiffuse/3.6,
	1                               fd,DT0)
           
CNB Note that need dTcorr ( %/degC : % increase for Deg C increase (typically negative) 
C                  and Tref (C) for this.
          dTcorr = Par(22) ! normally, this is the a_1 King parameter (a negative number normally)
          Tref   = Par(23) ! normally, this is the a_2 King Parameter
          Pmp = MAX(0.00, Pmp + Pmp*(dTcorr/100)*(Tcell-Tref))
        else
           tback = tamb
           tcell = tamb
        Endif
CNB In the case of single-pt efficiency the default operating voltage given to the model is used for all voltage outputs        
        Vmp = OperatingVoltage

C Other output values are set to zero and should not be used by SAM
         OUT(1)=OperatingVoltage       
         OUT(2)=0       
         OUT(3)=0       
 	   OUT(4)=EffMax  
         OUT(5)=Vmp       
         OUT(6)=0       
         OUT(7)=Pmp     
	   OUT(8)=EffMax  
         OUT(9)=0       
         OUT(10)=0      
         OUT(11)=Tcell      
         OUT(12)=Tback      
         OUT(13)=0     
         OUT(14)=0     
         OUT(15)=0    
	   OUT(16)=0      
	   OUT(17)=0    
	   OUT(18)=0     
	   OUT(19)=0     
         OUT(20)=1.0
         OUT(21)=0     
         OUT(22)=0      
         OUT(23)=0     
         OUT(24)=0  
         OUT(25)=0             
         OUT(26)=0               
         OUT(27)=0            
         OUT(28)=0             
         OUT(29)=0              
	   OUT(30)=0     
	   OUT(31)=0  
	   OUT(32)=0    
	   OUT(33)=0    
	   OUT(34)=0   
	   OUT(35)=0   
	   OUT(36)=0 
      
      Return 1
      Endif


CNB 12-23-04 CONCENTRATING PV: Single-efficiency model option put in to allow simple analysis with SAM
      IF ((a_0.LT.0).AND.(a_3.LT.0)) Then
CNB    a_0 is normally the negative fractional efficiency of the PV panel

C   Linearly interpolate efficiency from radiation level         
         TotalRad = (IcBeam)/3.6
         EffMax = EffInterpolate( TotalRad, MpeRadArr, MpeEffArr )
         
C  Pmp is the only real output available from this simple model
CNB This is only for a single module        
         Pmp = (EffMax*IcBeam*AColl)*NumStrings*NumSeries

CNB 2-20-10 added temperature correction to the CPV calcuations. These temp corrections exactly match the simple-efficiency and Sandia module
C           calculations as well. Diffuse radiation is set to 0 however. 

        If (Pmp > 0.001) THEN
CNB 7-25-05 Adding code to correct for temperature correction
CNB use King method to calculate cell temperature
CNB Note that need appropriate fd, a and b for this call
C Calculate back-of-module temperature:
          Tback = SandiaModuleTemperature(IcBeam/3.6,IcDiffuse/3.6,
	1                                WindSpeed,Tamb,fd,a,b)

CNB Note that need fd and DT0 for this call
C Calculate cell temperature:
          Tcell = SandiaTcellFromTmodule(Tback,IcBeam/3.6,
     !             IcDiffuse/3.6,fd,DT0)
           
CNB Note that need dTcorr ( %/degC : % increase for Deg C increase (typically negative) 
C                  and Tref (C) for this.
          dTcorr = Par(22) ! normally, this is the a_1 King parameter (a negative number normally)
          Tref   = Par(23) ! normally, this is the a_2 King Parameter
          Pmp = MAX(0.00, Pmp + Pmp*(dTcorr/100)*(Tcell-Tref))
        else
           tback = tamb
           tcell = tamb
        Endif
CNB In the case of CPV the default operating voltage given to the model is used for all voltage outputs        
        Vmp = OperatingVoltage

C Other output values are set to zero and should not be used by SAM
         OUT(1)=OperatingVoltage       
         OUT(2)=0       
         OUT(3)=0       
 	   OUT(4)=EffMax  
         OUT(5)=Vmp       
         OUT(6)=0       
         OUT(7)=Pmp     
	   OUT(8)=EffMax  
         OUT(9)=0       
         OUT(10)=0      
         OUT(11)=Tcell      
         OUT(12)=Tback      
         OUT(13)=0     
         OUT(14)=0     
         OUT(15)=0    
	   OUT(16)=0      
	   OUT(17)=0    
	   OUT(18)=0     
	   OUT(19)=0     
         OUT(20)=1.0
         OUT(21)=0     
         OUT(22)=0      
         OUT(23)=0     
         OUT(24)=0  
         OUT(25)=0             
         OUT(26)=0               
         OUT(27)=0            
         OUT(28)=0             
         OUT(29)=0              
	   OUT(30)=0     
	   OUT(31)=0  
	   OUT(32)=0    
	   OUT(33)=0    
	   OUT(34)=0   
	   OUT(35)=0   
	   OUT(36)=0 
      
      Return 1
      Endif




C Actual Sandia Model Starts Here

      IF ((IcBeam+IcDiffuse).GT.MinIrradiance) THEN

C Voltage across single module:
	  SingleModuleVoltage = OperatingVoltage/NumSeries

C Calculate back-of-module temperature:
        Tback = SandiaModuleTemperature(IcBeam/3.6,IcDiffuse/3.6,
	1                                WindSpeed,Tamb,fd,a,b)

C Calculate cell temperature:
        Tcell = SandiaTcellFromTmodule(Tback,IcBeam/3.6,IcDiffuse/3.6,
	1                               fd,DT0)

C Calculate Air Mass
        AMa = AbsoluteAirMass(ZenithAngle, Altitude)

C Calculate F1 function:
        F1 = SandiaF1(AMa,a_0,a_1,a_2,a_3,a_4)

C Calculate F2 function:
        F2 = SandiaF2(IncidenceAngle,b_0,b_1,b_2,b_3,b_4,b_5)

C Calculate short-circuit current:
        Isc = SandiaIsc(Tcell,Isc0,IcBeam/3.6,IcDiffuse/3.6,
	1                  F1,F2,fd,aIsc)

C Calculate effective irradiance:
        Ee = SandiaEffectiveIrradiance(Tcell,Isc,Isc0,aIsc)

C Calculate Imp:
        Imp = SandiaImp(Tcell,Ee,Imp0,aImp,c_0,c_1)

C Calculate Voc:
        Voc = SandiaVoc(Tcell,Ee,Voc0,NcellSer,DiodeFactor,
	1                BVoc0,mBVoc)

C Calculate Vmp:
        Vmp = SandiaVmp(Tcell,Ee,Vmp0,NcellSer,DiodeFactor,
	1                        BVmp0,mBVmp,c_2,c_3)

C Calculate Ix:
        Ix = SandiaIx(Tcell,Ee,Ix0,aIsc,aImp,c_4,c_5)

C Calculate Vx:
        Vx = Voc/2.0

C Calculate Ixx:
        Ixx = SandiaIxx(Tcell,Ee,Ixx0,aImp,c_6,c_7)

C Calculate Vxx:
        Vxx = 0.5*(Voc+Vmp)

C Calculate Pmp, single module:
        Pmp = Imp * Vmp * 3.6 !kJ/hr

C Calculate max efficiency:
        EffMax = Pmp/(IcBeam+IcDiffuse)/Acoll

C Attempt to fit single-diode model if Ic >= IrradianceThreshold;
C otherwise just use TRW equation.

C NB  7-20-05 SAM only uses the Maximum Power Point and so doesn't need to do the rest of this code.
C             Craig is trying to use some strange King values which crash during the following code
C             Therefore, this section is now skipped in SAM (probably still present in BE-Opt however)
      Goto 500

        IF (IcBeam+IcDiffuse.GE.IrradianceThreshold) THEN
      
C Calculate 5 parameters for single-diode model:

          MaxIter=200
          FitTolerance=1.0e-4
          Rs_Step=0.1
          z_Step=0.1
          Rs_InitialMatrixWidth=0.1
          z_InitialMatrixWidth=0.99
	    TRW_Ix=TRW_Current (Vx,Vmp,Voc,Isc,Imp)
	    TRW_Ixx=TRW_Current(Vxx,Vmp,Voc,Isc,Imp)
	    Rs_init=0.0
	    z_init=1.0

          CALL FindBestFitSingleDiode(Isc, Ix, Imp, Ixx,
	1                                Vx, Vmp, Vxx, Voc,
     1                                MaxIter, FitTolerance,
     1                                Rs_Step, z_Step,
     1                                Rs_InitialMatrixWidth,
     1                                z_InitialMatrixWidth,
     1                                TRW_Ix, TRW_Ixx,
     1                                Converged, MaxError,
     1                                Rs_init, z_init, TotIter,
     1                                Rs, z, Rsh, Io, IL)
        ELSE
	    MaxError = 999.
	    Converged =.FALSE.
C Set diode-model related variables to 999 for easy recognition:
          Rs_init   = 999.
		  z_init    = 999.
		  TotIter   = 999
          Rs        = 999.
		  z         = 999.
		  Rsh       = 999.
		  Io        = 999.
		  IL        = 999.
	  ENDIF

C If NOT Converged, calculate maximum fractional error using TRW equation
C and compare to MaxError returned from FindBestFitSingleDiode routine:
        IF (Converged.EQ..FALSE.) THEN
          Current_x=TRW_Current(Vx,Vmp,Voc,Isc,Imp)
          Current_mp=TRW_Current(Vmp,Vmp,Voc,Isc,Imp)
          Current_xx=TRW_Current(Vxx,Vmp,Voc,Isc,Imp)

	    MaxTRWerror=DMAX1(DABS(Current_x-Ix)/Ix,
     1                    DABS(Current_mp-Imp)/Imp,
     1                    DABS(Current_xx-Ixx)/Ixx)

	    IF (MaxTRWerror.GT.MaxError) THEN
            UseSingleDiodeModel=.TRUE.
	    ELSE
            UseSingleDiodeModel=.FALSE.
	      MaxError=MaxTRWerror
	    ENDIF
	  ELSE
          UseSingleDiodeModel=.TRUE.
	  ENDIF

C Calculate approximate operating current using TRW equation:
        I_TRW = TRW_Current(SingleModuleVoltage,Vmp,Voc,Isc,Imp)

	  IF (UseSingleDiodeModel.EQ..TRUE.) THEN
	    OperatingCurrent=SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,
	1                                SingleModuleVoltage,I_TRW)
	1                                *NumStrings
	  ELSE
C This is the next best guess if single-diode model didn't converge:
          OperatingCurrent = I_TRW*NumStrings

        ENDIF	

C Calculate Absolute Errors and Fractional Errors at Vx, Vmp, and Vxx:
C------------------------------------
	  IF (UseSingleDiodeModel.EQ..TRUE.) THEN
          Current_x = SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,Vx,
     1                        TRW_Current(Vx,Vmp,Voc,Isc,Imp))

          Current_mp = SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,Vmp,
     1                        TRW_Current(Vmp,Vmp,Voc,Isc,Imp))

          Current_xx = SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,Vxx,
     1                        TRW_Current(Vxx,Vmp,Voc,Isc,Imp))
	  ELSE
          Current_x =TRW_Current(Vx,Vmp,Voc,Isc,Imp)
          Current_mp =TRW_Current(Vmp,Vmp,Voc,Isc,Imp)
          Current_xx =TRW_Current(Vxx,Vmp,Voc,Isc,Imp)
	  ENDIF

C Absolute error in kJ/hr:
        AbsErrorVx=(Current_x-Ix)*Vx*3.6
        AbsErrorVmp=(Current_mp-Imp)*Vmp*3.6
        AbsErrorVxx=(Current_xx-Ixx)*Vxx*3.6

C Fractional error:
		FracErrorVx=(Current_x-Ix)/Ix
		FracErrorVmp=(Current_mp-Imp)/Imp
		FracErrorVxx=(Current_xx-Ixx)/Ixx
C------------------------------------


C Scale Absolute Errors to entire array:
        AbsErrorVx=AbsErrorVx*NumStrings*NumSeries
        AbsErrorVmp=AbsErrorVmp*NumStrings*NumSeries
        AbsErrorVxx=AbsErrorVxx*NumStrings*NumSeries

C Calculate operating power for whole array:
        OperatingPower = OperatingCurrent*OperatingVoltage*3.6

C Don't allow negative power (can happen if OperatingVoltage>Voc)
	  OperatingPower = DMAX1(OperatingPower,0.0)

C Calculate operating efficiency:
        Eff = OperatingPower/(IcBeam+IcDiffuse)
     &            /(Acoll*NumSeries*NumStrings)

CNB MOVED TO BOTTOM
CC Scale 5 points to NumStrings and NumSeries:
C        Pmp = Pmp*NumStrings*NumSeries
C	  Imp = Imp*NumStrings
C	  Vmp = Vmp*NumSeries
C	  Isc = Isc*NumStrings
C	  Voc = Voc*NumSeries
C	  Ix  = Ix*NumStrings
C	  Ixx = Ixx*NumStrings
C	  Vx  = Vx*NumSeries
C	  Vxx = Vxx*NumSeries
        

	ELSE ! Ibeam+Idiff < MaxIrradiance

	  OperatingCurrent=0.0
	  OperatingPower=0.0
	  Eff=0.0
	  Vmp=0.0
	  Imp=0.0
	  Pmp=0.0
	  EffMax=0.0
	  Isc=0.0
        Voc=0.0
        Tcell=Tamb
        Tback=Tamb
        AMa=999.0
        F1=0.0
        F2=0.0
	  Ix=0.0
	  Vx=0.0
	  Ixx=0.0
	  Vxx=0.0
	  Converged=.FALSE.
        MaxError=0.0
        Rs_init=0.0
        z_init=0.0
        TotIter=0.0
        Rs=0.0
        z=0.0
        Rsh=0.0
        Io=0.0
        IL=0.0
	  AbsErrorVx=0.0
	  AbsErrorVmp=0.0
	  AbsErrorVxx=0.0
	  FracErrorVx=0.0
	  FracErrorVmp=0.0
	  FracErrorVxx=0.0

	ENDIF !Ibeam+Idiff > MinIrradiance

500   Continue


CNB Copied down from above to get below the "GOTO bypass"
C Scale 5 points to NumStrings and NumSeries:
        Pmp = Pmp*NumStrings*NumSeries
	  Imp = Imp*NumStrings
	  Vmp = Vmp*NumSeries
	  Isc = Isc*NumStrings
	  Voc = Voc*NumSeries
	  Ix  = Ix*NumStrings
	  Ixx = Ixx*NumStrings
	  Vx  = Vx*NumSeries
	  Vxx = Vxx*NumSeries

      MinInverterCapacity = Imp0*Vmp0*InverterCapacityFactor*3.6

C SET OUTPUTS
      OUT(1)=OperatingVoltage !(Volts)
      OUT(2)=OperatingCurrent !(Amps)
      OUT(3)=OperatingPower   !(kJ/hr)
 	OUT(4)=Eff              !(unitless)
      OUT(5)=Vmp              !(Volts)
      OUT(6)=Imp              !(Amps)
      OUT(7)=Pmp              !(kJ/hr)
	OUT(8)=EffMax           !(unitless)
      OUT(9)=Isc              !(Amps)
      OUT(10)=Voc             !(Volts)
      OUT(11)=Tcell           !(deg C)
      OUT(12)=Tback           !(deg C)
      OUT(13)=AMa             !(unitless)
      OUT(14)=F1              !(unitless)
      OUT(15)=F2              !(unitless)
	OUT(16)=Ix              !(Amps)
	OUT(17)=Vx              !(Volts)
	OUT(18)=Ixx             !(Amps)
	OUT(19)=Vxx             !(Volts)
      IF (Converged.EQ..TRUE.) THEN
	  OUT(20)=1.0
	ELSE
	  OUT(20)=0.0           !(1=yes,0=no)
	ENDIF
      OUT(21)=MaxError        !(fraction)
      OUT(22)=Rs_init         !(ohms)
      OUT(23)=z_init          !(unitless)
      OUT(24)=DFLOAT(TotIter) !(unitless)
      OUT(25)=Rs              !(ohms)
      OUT(26)=z               !(unitless)
      OUT(27)=Rsh             !(ohms)
      OUT(28)=Io              !(amps)
      OUT(29)=IL              !(amps)
	OUT(30)=AbsErrorVx      !(kJ/hr)
	OUT(31)=AbsErrorVmp     !(kJ/hr)
	OUT(32)=AbsErrorVxx     !(kJ/hr)
	OUT(33)=FracErrorVx     !(kJ/hr)
	OUT(34)=FracErrorVmp    !(kJ/hr)
	OUT(35)=FracErrorVxx    !(kJ/hr)
	OUT(36)=MinInverterCapacity !(kJ/hr)
		
	
      RETURN 1
      END

C -------------------------------------------------------------------------------
C FUNCTIONS and SUBROUTINES
C -------------------------------------------------------------------------------
C -------------------------------------------------------------------------------
	FUNCTION TRW_Current(V,VmaxPow,Voc,Isc,ImaxPow)
C Returns current as a function of voltage using the TRW equation (A)
C VmaxPow = voltage at the maximum power point (V)
C ImaxPow = current at the maximum power point (A)
C Voc     = open-circuit voltage (V)
C Isc     = short-circuit current (A)

      IMPLICIT NONE

      DOUBLE PRECISION TRW_Current
	DOUBLE PRECISION V,VmaxPow,ImaxPow,Voc,Isc
	DOUBLE PRECISION C_1,C_2,Itrw
  
        IF ((Isc.GT.0.).AND.(Voc.GT.0.)) THEN
	    IF(ImaxPow.LT.Isc) THEN
            C_2 = (VmaxPow / Voc - 1.0)/ Log(1.0 - ImaxPow / Isc)
	    ELSE
	      C_2 = 0.0
	    ENDIF
          IF (C_2.GT.0.) THEN
            C_1 = (1.0 - ImaxPow / Isc) * Exp(-VmaxPow / C_2 / Voc)
      
            Itrw = Isc*(1.0 - C_1 * (Exp(V / C_2 / Voc) - 1.0))
          ELSE
            Itrw = 0.0
          ENDIF
        ELSE
          Itrw = 0.0
        ENDIF
        
	  IF(Itrw.LT.0.0) THEN
	    Itrw=0.0
	  ENDIF

	  TRW_Current=Itrw
      
      END FUNCTION
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaModuleTemperature(Ibc,Idc,Ws,Ta,fd,a,b)
C Returns back-of-module temperature, deg C
C Ibc = beam radiation on collector plane, W/m2
C Idc = Diffuse radiation on collector plane, W/m2
C Ws  = wind speed, m/s
C Ta  = ambient temperature, degC
C fd  = fraction of Idc used (empirical constant)
C a   = empirical constant
C b   = empirical constant

        IMPLICIT NONE

        DOUBLE PRECISION SandiaModuleTemperature
	  DOUBLE PRECISION Ibc,Idc,Ws,Ta,fd,a,b,E
  
        E = Ibc + fd * Idc
  
        SandiaModuleTemperature = E * Exp(a + b * Ws) + Ta
    
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaTcellFromTmodule(Tm,Ibc,Idc,fd,DT0)
C Returns cell temperature, deg C
C Tm  = module temperature (deg C)
C Ibc = beam radiation on collector plane, W/m2
C Idc = Diffuse radiation on collector plane, W/m2
C fd  = fraction of Idc used (empirical constant)
C DT0 = (Tc-Tm) at E=1000 W/m2 (empirical constant known as dTc), deg C

        IMPLICIT NONE
		
		DOUBLE PRECISION SandiaTcellFromTmodule
		DOUBLE PRECISION Tm,Ibc,Idc,fd,DT0
		DOUBLE PRECISION E
  
        E = Ibc + fd * Idc
    
        SandiaTcellFromTmodule = Tm + E / 1000.0 * DT0
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaCellTemperature(Ibc,Idc,Ws,Ta,fd,a,b,DT0)
C Returns cell temperature, deg C
C Ibc = beam radiation on collector plane, W/m2
C Idc = Diffuse radiation on collector plane, W/m2
C Ws  = wind speed, m/s
C Ta  = ambient temperature, degC
C fd  = fraction of Idc used (empirical constant)
C a   = empirical constant
C b   = empirical constant
C DT0 = (Tc-Tm) at E=1000 W/m2 (empirical constant known as dTc), deg C

        DOUBLE PRECISION SandiaCellTemperature
	  DOUBLE PRECISION Ibc,Idc,Ws,Ta,fd,a,b,DT0,E,Tm
  
        E = Ibc + fd * Idc
  
        Tm = E * Exp(a + b * Ws) + Ta
  
        SandiaCellTemperature = Tm + E / 1000.0 * DT0
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaEffectiveIrradiance(Tc,Isc,Isc0,aIsc)
C Returns "effective irradiance", used in calculation of Imp, Voc, Ix, Ixx
C Tc   = cell temperature
C Isc = short-circuit current under operating conditions (A)
C Isc0 = Isc at Tc=25 C, Ic=1000 W/m2 (A)
C aIsc = Isc temperature coefficient (degC^-1)
        DOUBLE PRECISION SandiaEffectiveIrradiance
	  DOUBLE PRECISION Tc,Isc,Isc0,aIsc

        SandiaEffectiveIrradiance = Isc / (1.0+aIsc*(Tc - 25.0))/Isc0
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function AbsoluteAirMass(SolZen, Altitude)
C Returns absolute air mass
C SolZen = solar zenith (deg)
C Altitude = site altitude (m)

        DOUBLE PRECISION AbsoluteAirMass
	  DOUBLE PRECISION SolZen, Altitude,AM
  
        IF (SolZen.LT.89.9) THEN
          AM = (Cos(SolZen * 0.01745329) + 0.5057 *
	1           (96.08 - SolZen)**-1.634)**-1.0
  
          AbsoluteAirMass = Exp(-0.0001184 * Altitude) * AM
        ELSE
          AbsoluteAirMass = 999.
        ENDIF

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaF1(AMa,a0,a1,a2,a3,a4)
C Returns the Sandia Air Mass function
C AMa = absolute air mass
C a0,a1,a2,a3,a4 = empirical constants, module-specIFic
  
        DOUBLE PRECISION SandiaF1
	  DOUBLE PRECISION AMa,a0,a1,a2,a3,a4,F1

        F1 = a0 + a1*AMa + a2*AMa**2.0 + a3*AMa**3.0 + a4*AMa**4.0
  
        IF (F1.GT.0.0) THEN
          SandiaF1 = F1
        ELSE
          SandiaF1 = 0.0
        ENDIF
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaF2(IncAng,b0,b1,b2,b3,b4,b5)
C Returns Sandia F2 function
C IncAng = incidence angle (deg)
C b0,b1,b2,b3,b4,b5 = empirical module-specIFic constants

        DOUBLE PRECISION SandiaF2
	  DOUBLE PRECISION IncAng,b0,b1,b2,b3,b4,b5,F2
  
        F2 = b0+b1*IncAng+b2*IncAng**2.0+b3*IncAng**3.0+b4*IncAng**4.0+
	1       b5*IncAng**5.0
  
        IF (F2.GT.0.0) THEN
          SandiaF2 = F2
        ELSE
          SandiaF2 = 0.0
        ENDIF
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaImp(Tc,Ee,Imp0,aImp,C0,C1)
C Returns current at maximum power point (A)
C Tc = cell temperature (degC)
C Ee = effective irradiance (W/m2)
C Imp0 = current at MPP at SRC (1000 W/m2, 25 C) (A)
C aImp = Imp temperature coefficient (degC^-1)
C c0,c1 = empirical module-specIFic constants

        DOUBLE PRECISION SandiaImp
	  DOUBLE PRECISION Tc,Ee,Imp0,aImp,C0,C1

        SandiaImp = Imp0*(C0*Ee+C1*Ee**2.0)*(1.0+aImp*(Tc-25.0))
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaIsc(Tc,Isc0,Ibc,Idc,F1,F2,fd,aIsc)
C Returns Short-Circuit Current
C Isc0 = Isc at Tc=25 C, Ic=1000 W/m2 (A)
C Ibc  = beam radiation on collector plane (W/m2)
C Idc  = Diffuse radiation on collector plane (W/m2)
C F1   = Sandia F1 function of air mass
C F2   = Sandia F2 function of incidence angle
C fd   = module-specIFic empirical constant
C aIsc = Isc temperature coefficient (1/degC)
C Tc   = cell temperature

        IMPLICIT NONE
        DOUBLE PRECISION SandiaIsc
	  DOUBLE PRECISION Tc
	  DOUBLE PRECISION Isc0
	  DOUBLE PRECISION Ibc
	  DOUBLE PRECISION Idc
	  DOUBLE PRECISION F1
	  DOUBLE PRECISION F2
	  DOUBLE PRECISION fd
	  DOUBLE PRECISION aIsc

CNB UPDATED 12-3-07 BASED ON INFO FROM UPDATED CODE FROM GREG BARKER
C        SandiaIsc=Isc0*((Ibc*F1*F2+fd*Idc)/1000.0)*(1.0+aIsc*(Tc-25.0))
        SandiaIsc=Isc0*F1*((Ibc*F2+fd*Idc)/1000.0)*(1.0+aIsc*(Tc-25.0))
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaIx(Tc,Ee,Ix0,aIsc,aImp,C4,C5)
C Returns current "Ix" at V=0.5*Voc (A)
C Tc = cell temperature (deg C)
C Ee = effective irradiance
C Ix0 = Ix at SRC (1000 W/m2, 25 C) (A)
C aIsc = Isc temp coefficient (/C)
C aImp = Imp temp coefficient (/C)
C c4,c5 = empirical module-specIFic constants

        DOUBLE PRECISION SandiaIx
	  DOUBLE PRECISION Tc,Ee,Ix0,aIsc,aImp,C4,C5

        SandiaIx=Ix0*(C4*Ee+C5*Ee**2.0)
	1               *(1.0+((aIsc+aImp)/2.0*(Tc-25.0)))

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaIxx(Tc,Ee,Ixx0,aImp,C6,C7)
C Returns current "Ix" at V=0.5*(Voc+Vmp) (A)
C Tc = cell temperature (deg C)
C Ee = effective irradiance
C Ixx0 = Ixx at SRC (1000 W/m2, 25 C) (A)
C aImp = Imp temp coefficient (/C)
C c4,c5 = empirical module-specIFic constants

        DOUBLE PRECISION SandiaIxx
	  DOUBLE PRECISION Tc,Ee,Ixx0,aImp,C6,C7

        SandiaIxx=Ixx0*(C6*Ee+C7*Ee**2.0)*(1.0+aImp*(Tc-25.0))

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaVmp(Tc,Ee,Vmp0,NcellSer,DiodeFactor,
	1                        BVmp0,mBVmp,C2,C3)
C Returns Voltage at Max. Power Point (V)
C Tc = cell temperature (deg C)
C Ee = effective irradiance
C Vmpo = Vmp at SRC (1000 W/m2, 25 C) (V)
C NcellSer = # cells in series
C DiodeFactor = module-specIFic empirical constant
C BVmp0 = Vmp temperature coefficient (V/C)
C mBVmp = change in BVmp with irradiance
C c2,c3 = empirical module-specIFic constants

        DOUBLE PRECISION SandiaVmp
	  DOUBLE PRECISION Tc,Ee,Vmp0,NcellSer,DiodeFactor,BVmp0,mBVmp,C2,C3
	  DOUBLE PRECISION dTc,BVmpEe
  
        IF (Ee.GT.0.0) THEN
          dTc=DiodeFactor*((1.38066E-23*(Tc+273.15))/1.60218E-19)
          BVmpEe = BVmp0 + mBVmp * (1.0 - Ee)

          SandiaVmp = Vmp0+C2*NcellSer*dTc*Log(Ee)+
	1                C3*NcellSer*(dTc*Log(Ee))**2+BVmpEe*(Tc-25.0)
        ELSE
          SandiaVmp = 0
        ENDIF

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SandiaVoc(Tc,Ee,Voc0,NcellSer,DiodeFactor,
	1                        BVoc0,mBVoc)
C Returns Open-Circuit Voltage (V)
C Tc = cell temperature (deg C)
C Ee = effective irradiance
C Voco = Voc at SRC (1000 W/m2, 25 C) (V)
C NcellSer = # cells in series
C DiodeFactor = module-specIFic empirical constant
C BVoc0 = Voc temperature coefficient (V/C)
C mBVoc = change in BVoc with irradiance

        DOUBLE PRECISION SandiaVoc
	  DOUBLE PRECISION Tc,Ee,Voc0,NcellSer,DiodeFactor,BVoc0,mBVoc
	  DOUBLE PRECISION dTc,BVocEe
  
        IF (Ee.GT.0.0) THEN
          dTc=DiodeFactor*((1.38066E-23*(Tc + 273.15))/1.60218E-19)
          BVocEe = BVoc0 + mBVoc * (1.0 - Ee)
  
          SandiaVoc=Voc0+NcellSer*dTc*Log(Ee)+BVocEe*(Tc-25.0)
        ELSE
          SandiaVoc = 0.0
        ENDIF
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function TRW_C1(VmaxPow,ImaxPow,Voc,Isc,C2)
C Returns C1 in the TRW equation for PVs
C VmaxPow = voltage at the maximum power point (V)
C ImaxPow = current at the maximum power point (A)
C Voc     = open-circuit voltage (V)
C Isc     = short-circuit current (A)
C C2      = C2 in the TRW equation for PVs

        DOUBLE PRECISION TRW_C1
	  DOUBLE PRECISION VmaxPow,ImaxPow,Voc,Isc,C2

        IF (((Isc.GT.0.).AND.(Voc.GT.0.)).AND.(C2.GT.0.)) THEN
          TRW_C1 = (1.0 - ImaxPow / Isc) * Exp(-VmaxPow / C2 / Voc)
        ELSE
          TRW_C1 = 0.0
        ENDIF
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function TRW_C2(VmaxPow,ImaxPow,Voc,Isc)
C Returns C2 in the TRW equation for PVs
C VmaxPow = voltage at the maximum power point (V)
C ImaxPow = current at the maximum power point (A)
C Voc     = open-circuit voltage (V)
C Isc     = short-circuit current (A)

        DOUBLE PRECISION TRW_C2
	  DOUBLE PRECISION VmaxPow,ImaxPow,Voc,Isc

        IF ((Isc.GT.0.).AND.(Voc.GT.0.)) THEN
          TRW_C2 = (VmaxPow / Voc - 1.0) / Log(1.0 - ImaxPow / Isc)
        ELSE
          TRW_C2 = 0.0
        ENDIF
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Subroutine FindBestFitSingleDiode(Isc, Ix, ImaxP, Ixx, Vx, VmaxP,
	1                                  Vxx, Voc, MaxIter, FitTolerance,
     1                                  Rs_Step, z_Step,
     1                                  Rs_InitialMatrixWidth,
     1                                  z_InitialMatrixWidth,
     1                                  TRW_Ix, TRW_Ixx,
     1                                  Converged, MaxError,
     1                                  Rs_init, z_init, TotIter,
     1                                  Rs, z, Rsh, Io, IL)
  
C INPUT:
C MaxIter               = maximum number of iterations within minimization routine
C FitTolerance          = RMS error at which fit is considered adequate
C Rs_Step               = step size for Rs when trying new initial value
C z_Step                = step size for z when trying new initial value
C Rs_InitialMatrixWidth = initial matrix of RMS values will be generated with values
C                         of Rs=Rs_init (+/-) Rs_InitialMatrixWidth
C z_InitialMatrixWidth  = initial matrix of RMS values will be generated with values
C                         of z=z_init (+/-) z_InitialMatrixWidth
C TRW_Ix                = Ix as calculated by the TRW equation
C TRW_Ixx               = Ixx as calculated by the TRW equation

C RETURNED:
C Converged             = boolean variable returned indicating whether or not the
C                         solution converged to the desired tolerance
C MaxError              = maximum error (fraction) using final values of Rs and z
C                         in single-diode equation to calculate Ix, Ixx
C Rs_init               = initial value of Rs fed to minimization routine (the value
C                         of Rs_init input as a parameter is the initial value of Rs
C                         used in the solution). Recommend Rs_init=0.0
C z_init                = initial value of z fed to minimization routine (the value
C                         of z_init input as a parameter is the initial value of z
C                         used in the solution). Recommend z_init=1.0
C TotIter               = total number of iterations made
C Rs                    = series resistance
C z                     = fitting parameter (the value of z input as a parameter
C                         is the initial value of z used in the solution). Recommend z=1.0
C Rsh                   = shunt resistance
C Io                    = diode current
C IL                    = light current
  
C Subroutine parameters:      
      DOUBLE PRECISION Isc, Ix, ImaxP, Ixx, Vx, VmaxP, Vxx, Voc
	INTEGER MaxIter
	DOUBLE PRECISION FitTolerance
      DOUBLE PRECISION Rs_Step, z_Step
      DOUBLE PRECISION Rs_InitialMatrixWidth, z_InitialMatrixWidth
      DOUBLE PRECISION TRW_Ix, TRW_Ixx
      LOGICAL Converged
	DOUBLE PRECISION MaxError
      DOUBLE PRECISION Rs_init, z_init
	INTEGER TotIter
      DOUBLE PRECISION Rs, z, Rsh, Io, IL

C Internal variables:
      INTEGER i,Niter
      DOUBLE PRECISION BestError,Best_RsInit,Best_zInit
      DOUBLE PRECISION Rs_last
      DOUBLE PRECISION z_last
      DOUBLE PRECISION OriginalRs_init
	INTEGER ErrorCode
  
      DOUBLE PRECISION Y(3)
      DOUBLE PRECISION FracError(2)
      DOUBLE PRECISION InitVal(2)
      DOUBLE PRECISION Bounds(2,2)
      DOUBLE PRECISION P(3,2)
      DOUBLE PRECISION OtherParams(10)
  
      DOUBLE PRECISION Max_z
	DOUBLE PRECISION Max_Rs
      DOUBLE PRECISION MinimizeTolerance
    
        Max_z  = 5.0
	  Max_Rs = 3.0
        MinimizeTolerance = 0.0001
	  
	  TotIter = 0
  
        OriginalRs_init = Rs_init
  
        z_init = z_init - z_Step
        InitVal(2) = z_init
        Best_zInit = z_init
        z_last = 0.0
        Best_RsInit = OriginalRs_init
        
        MaxError = 100.0
        BestError = 100.0
  
        Converged = .TRUE.
  
        DO WHILE ((DABS(MaxError).GT.FitTolerance).AND.
     &           (z_last.LE.Max_z))

C Increment initial value of z:
          z_init = z_init + z_Step
      
C give up; go back to the best solution:
          IF ((Rs_init.GE.Max_Rs).AND.(z_init.GE.Max_z)) THEN
            z_init = Best_zInit
C set z_last to number such that WHILE loop will not continue:
            z_last = Max_z + 1.0
            Converged = .FALSE.
          ELSE
            z_last = z_init
            Rs_init = OriginalRs_init - Rs_Step
            Rs_last = Rs_init
          ENDIF
  
          DO WHILE ((DABS(MaxError).GT.FitTolerance)
	1               .AND.(Rs_last.LE.Max_Rs))

C Increment initial value of Rs:
            Rs_init = Rs_init + Rs_Step
    
C give up; go back to the best solution:
            IF ((Rs_init.GE.Max_Rs).AND.(z_init.GE.Max_z)) THEN
              Rs_init = Best_RsInit
              z_init = Best_zInit
C set Rs_last, z_last to numbers such that WHILE loop will not continue:
              Rs_last = Max_Rs + 0.1
              z_last = Max_z + 1.0
              Converged = .FALSE.
            ELSE
              Rs_last = Rs_init
            ENDIF
    
            InitVal(1) = Rs_init
            InitVal(2) = z_init
C Don't allow z=0 to avoid divide-by-zero in single-diode current equation:
            IF (InitVal(2).EQ.0.0) THEN
              InitVal(2) = 1E-20
            ENDIF
      
C ---------------------------------------------------------------
C Minimize using AMOEBA minimization routine

            Bounds(1, 1) = -10.0 !minimum value of Rs
            Bounds(1, 2) = 10.0  !maximum value of Rs
            Bounds(2, 1) = 0.1   !minimum value of z
            Bounds(2, 2) = 10.0  !maximum value of z
  
C Fill in initial-value matrix:
            P(1, 1) = InitVal(1)
            P(1, 2) = InitVal(2)
  
            P(2, 1) = P(1, 1) + Rs_InitialMatrixWidth
            P(2, 2) = P(1, 2) + z_InitialMatrixWidth
            P(3, 1) = P(1, 1) - Rs_InitialMatrixWidth
            P(3, 2) = P(1, 2) - z_InitialMatrixWidth

            DO i = 1,3
              Rsh = SingleDiodeShuntResistance(P(i, 1), P(i, 2),
	1	                   	                   Isc, Voc, VmaxP, ImaxP)
              Io = SingleDiodeDiodeCurrent(P(i, 1), P(i, 2),
	1				                             Isc, Voc, Rsh)
              IL = SingleDiodeLightCurrent(P(i, 2), Voc, Rsh, Io, Isc)
              Y(i) = RMSerr(P(i, 1), P(i, 2), Rsh, Io, IL,
	1                      TRW_Ix, TRW_Ixx,Ix,Ixx,Vx,Vxx)
            END DO
       
            OtherParams(1) = Vx
            OtherParams(2) = TRW_Ix
            OtherParams(3) = Vxx
            OtherParams(4) = TRW_Ixx
            OtherParams(5) = Isc
            OtherParams(6) = Voc
            OtherParams(7) = ImaxP
            OtherParams(8) = VmaxP
            OtherParams(9) = Ix
            OtherParams(10) = Ixx
       
C Call AMOEBA routine to find values of Rs and z which result in minimum RMS error:
            CALL AMOEBA (P,Y,Bounds,OtherParams,2,
	1                  MinimizeTolerance, MaxIter, Niter, ErrorCode)
    
            TotIter = TotIter + Niter
    
            Rs = P(2, 1)
            z = P(2, 2)
C Recalculate other three parameters:
            Rsh = SingleDiodeShuntResistance(Rs, z, Isc, Voc,
	1                                       VmaxP, ImaxP)
            Io = SingleDiodeDiodeCurrent(Rs, z, Isc, Voc, Rsh)
            IL = SingleDiodeLightCurrent(z, Voc, Rsh, Io, Isc)
C ---------------------------------------------------------------
  
C Fraction error for each of the two (x and xx) points:
            FracError(1) = 1.0 - SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,
	1                                             Vx,TRW_Ix) / Ix
            FracError(2) = 1.0 - SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,
	1                                             Vxx, TRW_Ixx) / Ixx
  
            MaxError = 0.0
            DO i = 1,2
              IF (DABS(FracError(i)).GT.MaxError) THEN
                MaxError = DABS(FracError(i))
              ENDIF
            END DO
  
            IF (MaxError.LT.BestError) THEN
              BestError = MaxError
              Best_RsInit = Rs_init
              Best_zInit = z_init
            ENDIF
  
          END DO
        END DO


  
      End Subroutine

C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      SUBROUTINE AMOEBA(P,Y,Bounds,OtherParams,NP,FTOL,ItMax,ITER,Ecode)

C Slightly altered from Numerical Recipes book.

C Multidimensional minimization of the function USERFUNC(x()) where x is an
C NDIM-dimensional vector, by the downhill simplex method of Nelder and Mead.
C Input is a matrix P whose NDIM+1 rows are NDIM-dimensional vectors which are the
C vertices of the starting simplex (Logical dimensions of P are P(NDIM+1,NDIM);
C physical dimensions are input as P(NP+1,NP) ). Also input is the vector Y of length
C NDIM+1, whose components must be pre-initialized to the values of FunctionToMinimize
C evaluated at the NDIM+1 vertices (rows) of P; and FTOL the fractional convergence tolerance
C to be achieved in the function value (n.b.!). On output, P and Y will have been
C reset to NDIM+1 new points all within FTOL of a minimum function value, and ITER
C gives the number of iterations taken. Ecode is returned as zero IF all is well;
C Ecode=1 means ITMAX was exceeded.
  
C For passing to FunctionToMinimize function:
C OtherParams(1)= Vx
C OtherParams(2)= TRW_Ix
C OtherParams(3)= Vxx
C OtherParams(4)= TRW_Ixx
C OtherParams(5)= Isc
C OtherParams(6)= Voc
C OtherParams(7)= ImaxP
C OtherParams(8)= VmaxP
C OtherParams(9)= Ix
C OtherParams(10)= Ixx

        PARAMETER (NMAX=20,ALPHA=1.0,BETA=0.5,GAMMA=2.0)
  
        DOUBLE PRECISION P
        DIMENSION P(NP+1,NP)

        DOUBLE PRECISION Y
        DIMENSION Y(NP+1)

	  DOUBLE PRECISION PR
        DIMENSION PR(NMAX)

        DOUBLE PRECISION PRR
        DIMENSION PRR(NMAX)

        DOUBLE PRECISION PBAR
        DIMENSION PBAR(NMAX)

        DOUBLE PRECISION Bounds
        DIMENSION Bounds(NP,NP)
        
        DOUBLE PRECISION OtherParams
        DIMENSION OtherParams(10)

	  INTEGER NP
        DOUBLE PRECISION FTOL
	  INTEGER ITMAX
	  INTEGER ITER
	  INTEGER Ecode
        
	  DOUBLE PRECISION Denom
        INTEGER NDIM
	  INTEGER MPTS
	  INTEGER i,j,k
	  INTEGER ILO,IHI,INHI
	  DOUBLE PRECISION RTOL

        NDIM = NP
  
        Ecode = -1
  
        MPTS = NDIM + 1
        ITER = 0
  
        DO WHILE (Ecode.LT.0)
C ---------------------------------------------------------------------
C Don't allow variables to go out of bounds
          DO i = 1,NP
            DO j = 1,NP + 1
              IF (P(j, i).LT.Bounds(i, 1)) THEN
                P(j, i) = Bounds(i, 1)
              ENDIF
              IF (P(j, i).GT.Bounds(i, 2)) THEN
                P(j, i) = Bounds(i, 2)
              ENDIF
            END DO
          END DO
C ---------------------------------------------------------------------
          ILO = 1
          IF (Y(1).GT.Y(2)) THEN
            IHI = 1
            INHI = 2
          ELSE
            IHI = 2
            INHI = 1
          ENDIF
      
          DO i = 1,MPTS
            IF (Y(i).LT.Y(ILO)) THEN
	        ILO = i
	      ENDIF
      
	      IF (Y(i).GT.Y(IHI)) THEN
              INHI = IHI
              IHI = i
            ELSE
              IF (Y(i).GT.Y(INHI)) THEN
                IF (i.NE.IHI) THEN
	            INHI = i
	          ENDIF
              ENDIF
            ENDIF
          END DO
      
          Denom = (DABS(Y(IHI)) + DABS(Y(ILO)))
          IF (Denom.EQ.0.0) THEN Denom = 1E-20
    
          RTOL = 2.0 * DABS(Y(IHI) - Y(ILO)) / Denom
      
          IF (RTOL.LT.FTOL) THEN
            Ecode = 0
            RETURN
          ENDIF
      
          IF (ITER.EQ.ItMax) THEN
            Ecode = 1
          ENDIF
      
          ITER = ITER + 1
          DO j = 1,NDIM
            PBAR(j) = 0.0
          END DO
      
          DO i = 1,MPTS
            IF (i.NE.IHI) THEN
              DO j = 1,NDIM
                PBAR(j) = PBAR(j) + P(i, j)
              END DO
            ENDIF
          END DO
      
          DO j = 1,NDIM
            PBAR(j) = PBAR(j) / NDIM
            PR(j) = (1.0 + ALPHA) * PBAR(j) - ALPHA * P(IHI, j)
          END DO
      
C ---------------------------------------------------------------------
C Don't allow variables to go out of bounds
          DO k = 1,NDIM
            IF (PR(k).LT.Bounds(k, 1)) THEN
              PR(k) = Bounds(k, 1)
            ENDIF
            IF (PR(k).GT.Bounds(k, 2)) THEN
              PR(k) = Bounds(k, 2)
            ENDIF
          END DO
C ---------------------------------------------------------------------
          YPR = FunctionToMinimize(PR,OtherParams)
     
          IF (YPR.LE.Y(ILO)) THEN
            DO j = 1,NDIM
              PRR(j) = GAMMA * PR(j) + (1.0 - GAMMA) * PBAR(j)
            END DO
C ---------------------------------------------------------------------
C Don't allow variables to go out of bounds
          DO k = 1,NDIM
            IF (PRR(k).LT.Bounds(k, 1)) THEN
              PRR(k) = Bounds(k, 1)
            ENDIF
            IF (PRR(k).GT.Bounds(k, 2)) THEN
              PRR(k) = Bounds(k, 2)
            ENDIF
          END DO
C ---------------------------------------------------------------------
            YPRR = FunctionToMinimize(PRR,OtherParams)
            IF (YPRR.LT.Y(ILO)) THEN
              DO j = 1,NDIM
                P(IHI, j) = PRR(j)
              END DO
              Y(IHI) = YPRR
            ELSE
              DO j = 1,NDIM
                P(IHI, j) = PR(j)
              END DO
              Y(IHI) = YPR
            ENDIF
          ELSE
            IF (YPR.GE.Y(INHI)) THEN
              IF (YPR.LT.Y(IHI)) THEN
                DO j = 1,NDIM
                  P(IHI, j) = PR(j)
                END DO
                Y(IHI) = YPR
              ENDIF
              
              DO j = 1,NDIM
                PRR(j) = BETA * P(IHI, j) + (1.0 - BETA) * PBAR(j)
              END DO
C ---------------------------------------------------------------------
C Don't allow variables to go out of bounds
          DO k = 1,NDIM
            IF (PRR(k).LT.Bounds(k, 1)) THEN
              PRR(k) = Bounds(k, 1)
            ENDIF
            IF (PRR(k).GT.Bounds(k, 2)) THEN
              PRR(k) = Bounds(k, 2)
            ENDIF
          END DO
C ---------------------------------------------------------------------
               YPRR = FunctionToMinimize(PRR,OtherParams)
              IF (YPRR.LT.Y(IHI)) THEN
                DO j = 1,NDIM
                  P(IHI, j) = PRR(j)
                END DO
                Y(IHI) = YPRR
              ELSE
                DO i = 1,MPTS
                  IF (i.NE.ILO) THEN
                    DO j = 1,NDIM
                      PR(j) = 0.5 * (P(i, j) + P(ILO, j))
                      P(i, j) = PR(j)
                    END DO
C ---------------------------------------------------------------------
C Don't allow variables to go out of bounds
          DO k = 1,NDIM
            IF (PR(k).LT.Bounds(k, 1)) THEN
              PR(k) = Bounds(k, 1)
            ENDIF
            IF (PR(k).GT.Bounds(k, 2)) THEN
              PR(k) = Bounds(k, 2)
            ENDIF
          END DO
C ---------------------------------------------------------------------
                    Y(i) = FunctionToMinimize(PR,OtherParams)
                  ENDIF
                END DO
              ENDIF
            ELSE
              DO j = 1,NDIM
                P(IHI, j) = PR(j)
              END DO
              Y(IHI) = YPR
            ENDIF
    
          ENDIF

        END DO

      END SUBROUTINE
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function FunctionToMinimize(X,Constants)

C Returns RMS error of fit to five points on IV curve

        DOUBLE PRECISION FunctionToMinimize
	  DOUBLE PRECISION Rs
        DOUBLE PRECISION z
        DOUBLE PRECISION Vx
        DOUBLE PRECISION TRW_Ix
        DOUBLE PRECISION Vxx
        DOUBLE PRECISION TRW_Ixx
        DOUBLE PRECISION Isc
        DOUBLE PRECISION Voc
        DOUBLE PRECISION ImaxP
        DOUBLE PRECISION VmaxP
        DOUBLE PRECISION Rsh
        DOUBLE PRECISION Io
        DOUBLE PRECISION IL
        DOUBLE PRECISION Ix_calc
        DOUBLE PRECISION Ixx_calc
        DOUBLE PRECISION Ix
	  DOUBLE PRECISION Ixx

        DOUBLE PRECISION AbsoluteError(2)
	  DOUBLE PRECISION MaxError
        DOUBLE PRECISION Weight(2)
        DOUBLE PRECISION SumOfWeights
        DOUBLE PRECISION RMSerr
        INTEGER i

	  DOUBLE PRECISION X(2)
	  DOUBLE PRECISION Constants(10)
    
        Rs = X(1)
        z = X(2)
        Vx = Constants(1)
        TRW_Ix = Constants(2)
        Vxx = Constants(3)
        TRW_Ixx = Constants(4)
        Isc = Constants(5)
        Voc = Constants(6)
        ImaxP = Constants(7)
        VmaxP = Constants(8)
        Ix = Constants(9)
        Ixx = Constants(10)

        Rsh = SingleDiodeShuntResistance(Rs, z, Isc, Voc, VmaxP, ImaxP)
        Io = SingleDiodeDiodeCurrent(Rs, z, Isc, Voc, Rsh)
        IL = SingleDiodeLightCurrent(z, Voc, Rsh, Io, Isc)
        Ix_calc = SingleDiodeCurrentTRW(IL, Io, Rs, z, Rsh, Vx, TRW_Ix)
        Ixx_calc = SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,Vxx,TRW_Ixx)
  
        AbsoluteError(1) = DABS(Ix - Ix_calc) / Ix
        AbsoluteError(2) = DABS(Ixx - Ixx_calc) / Ixx
  
        MaxError = 0.0
        DO i = 1,2
          IF (AbsoluteError(i).GT.MaxError) THEN
	      MaxError = AbsoluteError(i)
	    ENDIF
        END DO
  
        Weight(1) = 1.0
        Weight(2) = 1.0
  
        RMSerr = 0.0
        SumOfWeights = 0.0
        DO i = 1,2
          RMSerr = RMSerr + AbsoluteError(i)**2 * Weight(i)
          SumOfWeights = SumOfWeights + Weight(i)
        END DO
  
        FunctionToMinimize = (RMSerr / SumOfWeights)**0.5

      End Function

C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SingleDiodeShuntResistance(Rs,z,Isc,Voc,VmaxP,ImaxP)
  
        DOUBLE PRECISION SingleDiodeShuntResistance
	  DOUBLE PRECISION Rs,z,Isc,Voc,VmaxP,ImaxP
		
	  DOUBLE PRECISION Term1
        DOUBLE PRECISION Term2
        DOUBLE PRECISION Term3
        DOUBLE PRECISION Numerator
        DOUBLE PRECISION Denominator
  
        IF(z.EQ.0.0) THEN
          z = 1E-20
        ENDIF
  
        IF((Voc / z).GT.50.0) THEN
          Term1 = 50.0
        ELSE
          Term1 = Voc / z
        ENDIF
  
        IF(((VmaxP + ImaxP * Rs) / z).GT.50.0) THEN
          Term2 = 50.0
        ELSE
          Term2 = (VmaxP + ImaxP * Rs) / z
        ENDIF
  
        IF((Isc * Rs / z).GT.50.0) THEN
          Term3 = 50.0
        ELSE
          Term3 = Isc * Rs / z
        ENDIF
  
        Numerator = ((Isc * Rs - Voc) * (Exp(Term1) - Exp(Term2)) + 
	1        (Voc - VmaxP - ImaxP * Rs) * (Exp(Term1) - Exp(Term3)))

        Denominator = (ImaxP * (Exp(Term1) - Exp(Term3)) + 
	1         Isc * (Exp(Term2) - Exp(Term1)))
        
		IF(DABS(Denominator).LT.1E-30) THEN
          SingleDiodeShuntResistance = 1E+20
        ELSE
          SingleDiodeShuntResistance = Numerator / Denominator
        ENDIF
  
      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SingleDiodeDiodeCurrent(Rs, z, Isc, Voc, Rsh)
  
        DOUBLE PRECISION SingleDiodeDiodeCurrent
	  DOUBLE PRECISION Rs,z,Isc,Voc,Rsh
        DOUBLE PRECISION Term1,Term2

	  Term1=DMIN1(Voc/z,50.0)
	  Term2=DMIN1(Isc*Rs/z,50.0)
	  Rsh=DMAX1(Rsh,1e-20)

        SingleDiodeDiodeCurrent = (Isc * Rsh + Isc * Rs - Voc) / Rsh / 
	1                          	(Exp(Term1) - Exp(Term2))
  

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SingleDiodeLightCurrent(z, Voc, Rsh, Io, Isc)

	  DOUBLE PRECISION SingleDiodeLightCurrent
	  DOUBLE PRECISION z,Voc,Rsh,Io,Isc
        DOUBLE PRECISION Term1

	  IF (z.eq.0.0) THEN
	    z=1e-20
	  ENDIF
	  Term1=DMIN1(Voc/z,50.0)
	  Rsh=DMAX1(Rsh,1e-20)
  
  
        SingleDiodeLightCurrent = Voc / Rsh + Io * (Exp(Term1) - 1)
  

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function SingleDiodeCurrentTRW(IL, Io, Rs, z, Rsh, V, I_TRW)
C Returns current as a function of voltage, the five single-diode model parameters,
C and the TRW equation value of current.
    
        DOUBLE PRECISION SingleDiodeCurrentTRW
	  DOUBLE PRECISION IL,Io,Rs,z,Rsh,V,I_TRW
	  DOUBLE PRECISION Term1
    
        Term1 = DMIN1(((V + I_TRW * Rs) / z),50.0)

        SingleDiodeCurrentTRW = IL-Io*(Exp(Term1)-1)-(V+I_TRW*Rs)/Rsh

      End Function
C -------------------------------------------------------------------------------

C -------------------------------------------------------------------------------
      Function RMSerr(Rs,z,Rsh,Io,IL,TRW_Ix,TRW_Ixx,Ix,Ixx,Vx,Vxx)

C Returns RMS error of fit to five points on IV curve
        
	  DOUBLE PRECISION RMSerr
	  DOUBLE PRECISION Rs, z, Rsh, Io, IL, TRW_Ix, TRW_Ixx 
	  DOUBLE PRECISION Ix, Ixx, Vx, Vxx

        DOUBLE PRECISION Ix_calc
        DOUBLE PRECISION Ixx_calc
        DOUBLE PRECISION AbsoluteError(2)
        DOUBLE PRECISION Weight(2)
        DOUBLE PRECISION SumOfWeights
	  DOUBLE PRECISION MaxError
        INTEGER i
    
        Ix_calc = SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,Vx, TRW_Ix)
        Ixx_calc = SingleDiodeCurrentTRW(IL,Io,Rs,z,Rsh,Vxx,TRW_Ixx)
  
        AbsoluteError(1) = DABS(Ix - Ix_calc) / Ix
        AbsoluteError(2) = DABS(Ixx - Ixx_calc) / Ixx
  
        MaxError = 0.0
        DO i = 1,2
          IF(AbsoluteError(i).GT.MaxError) THEN
	      MaxError = AbsoluteError(i)
	    ENDIF
        END DO
  
        Weight(1) = 1.0
        Weight(2) = 1.0
  
        RMSerr = 0.0
        SumOfWeights = 0.0
        DO i = 1,2
          RMSerr = RMSerr + AbsoluteError(i)**2 * Weight(i)
          SumOfWeights = SumOfWeights + Weight(i)
        END DO
  
        RMSerr = (RMSerr / SumOfWeights)**0.5

      End Function      
      
C****************************************************************
      DOUBLE PRECISION FUNCTION EffInterpolate( RADLEV, RADARR, EFFARR )

C     added by Aron Dobos 25 sept 09 SAM development
C     linearly interpolate efficiency from 5 rad/eff pairs

      IMPLICIT NONE
      DOUBLE PRECISION RADLEV, RADARR(5), EFFARR(5),wx
      INTEGER i,i1
       
       IF (RADLEV.le.RADARR(1)) THEN
            EffInterpolate = EFFARR(1)
       ELSE IF (RADLEV.ge.RADARR(5)) THEN
            EffInterpolate = EFFARR(5)
       ELSE
       
C    Interpolate linearly
              DO i=2,5
                IF (RADLEV.lt.RADARR(i)) EXIT
              ENDDO
              i1=i-1
              wx=(RADLEV-RADARR(i1))/(RADARR(i)-RADARR(i1))
              EffInterpolate=(1-wx)*EFFARR(i1)+wx*EFFARR(i) 
       ENDIF
      
      END FUNCTION