Problems scripting the "having centroid in" spatial query....

[AlistairTowers]

I have tried to bring together some previous coding examples gleaned from other posts to emulate the Polygon to Polygon "having centroid in" spatial query. I'm struggling to get this to work, has anyone got some comments/suggestions on my code (attached below)?

FYI: Layer 0 is the Trade Area layer (built from amalgamated collection district boundaries) and Layer 1 is the census collection boundaries themselves. I want to build a list of the CD's within the TA so that I can then extract and summarise associated census information for each trade area sector.

Any assistance greatly appreciated.

Alistair

Public Sub CDCentroids_WithinTA()

     Dim Doc As IMxDocument
     Dim pMap As IMap
     Dim pTargetLayer As IFeatureLayer
     Dim pSourceLayer As IFeatureLayer
     Dim pSGeometryCollection As IGeometryCollection
     Dim pTGeometryCollection As IGeometryCollection
     Dim pSpatialFilter As ISpatialFilter
     Dim pSFeatureSelection As IFeatureSelection
     
     
     Set Doc = ThisDocument
     Set pMap = Doc.FocusMap
     Set pTargetLayer = pMap.Layer(0)
     Set pSourceLayer = pMap.Layer(1)
     
     
     Dim pTFC As IFeatureClass
     Set pTFC = pTargetLayer.FeatureClass
    
     Dim pSFC As IFeatureClass
     Set pSFC = pSourceLayer.FeatureClass
     
     Dim pFCursor As IFeatureCursor
     Set pFCursor = pTFC.Update(Nothing, False)
    
    
     Dim pSCursor As IFeatureCursor
     Set pSCursor = pSFC.Update(Nothing, False)
     
     'pFeatureSelection.SelectionSet.Search Nothing, False, pFCursor
    
     Dim pFeature As IFeature
     Set pFeature = pFCursor.NextFeature
     
     Dim pSFeature As IFeature
     Set pSFeature = pSCursor.NextFeature
     
     Set pSFeatureSelection = pSourceLayer
     Set pSGeometryCollection = New GeometryBag
     
     Dim pSArea As IArea
     Dim pPoint As IPoint
    
     
     Do While Not pSFeature Is Nothing
            
            'MsgBox "Sector: " & pSFeature.Value(pSFeature.Fields.FindField("CD_CODE"))
            Set pSArea = pSFeature.Shape
            Set pPoint = pSArea.Centroid
            'pSGeometryCollection.AddGeometry pSFeature.Shape
            pSGeometryCollection.AddGeometry pPoint
            Set pSFeature = pSCursor.NextFeature
     Loop
        
     
     
     
     
     
     
     Dim pRelOp As IRelationalOperator
     Set pRelOp = pTGeometryCollection
     
     'MsgBox "Sector: " & pFeature.Value(pFeature.Fields.FindField("TA_NAME"))
        
        Do While Not pFeature Is Nothing
            Set pTGeometryCollection = New GeometryBag
            pTGeometryCollection.AddGeometry pFeature.Shape
            Set pRelOp = pTGeometryCollection
            'MsgBox "Sector: " & pFeature.Value(pFeature.Fields.FindField("TA_NAME"))
       
            If (pRelOp.Contains(pSGeometryCollection)) Then
                pSFeatureSelection.Add pFeature
            End If
            
            Set pFeature = pFCursor.NextFeature
        Loop
    
    
          
     Doc.ActivatedView.Refresh
     
       
End Sub